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ANISOTROPY OF THE SPECTRUM OF TURBULENCE 
AT SMALL WAVE-NUMBERS 
3y G. K. BATCHELOR and R. W. STEWART 
(Cavendish Laboratory, Cambridge) 
Received 29 March 1949] 


SUMMARY 

In a previous paper it has been shown theoretically that the tendency to isotropic 
turbulence does not exist for the large-scale components of homogeneous turbulence. 
This paper provides experimental evidence that the large-scale components of turbu- 
lence behind a grid in a uniform stream are anisotropic, even though the total 
turbulent energy may be distributed with approximate spherical symmetry. The 
evidence consists (a) of a comparison between two different velocity correlations for 
large spatial intervals, in the initial period of decay, and (b) of measurements of the 
directional distribution of turbulent energy in the final period of decay when all 
except the large-scale components of the motion have been greatly reduced by 
viscous decay. 


1. Introduction 

It is a matter of common observation that there exists a strong tendency 
to isotropy in all fields of turbulence. When there are no opposing influences 
such as a non-uniform distribution of mean velocity, isotropic turbulence 
is quickly set up. This process has been used for many years as a means 
of producing isotropic turbulence for experimental work in a wind-tunnel. 
A grid of regularly spaced bars is placed in a uniform stream of air, and at 
a comparatively short distance downstream from the grid the series of 
wakes from the various bars has been transformed into isotropic turbu- 
lence superimposed on a mean velocity distribution which is once again 
uniform. 

This tendency to spherical symmetry of the turbulence is essentially a 
consequence of the action of the pressure terms in the Navier-Stokes 
equations. The static pressure at any point is a non-directional quantity 
and provides a means of converting energy associated with one component 
of velocity into energy associated with another component. 

Now in another paper (1) one of us has investigated the way in which 
the tendency to isotropy, for a particular Fourier component of the motion, 
varies with the wave-number of the component in a field of homogeneous 


turbulence. It was shown there that the energy spectrum, defined as the 


density in three-dimensional wave-number space of contributions to the 


energy of the turbulence, is of order k? (k magnitude of wave-number 
vector) when / is small, and that the rate of change with time of this 
energy spectrum is of higher order in k. When k is small, all three dynamical 
[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 0)] £ 
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effects represented in the Navier-Stokes equations—viscous decay, the 
tendency to isotropy owing to the pressure term, and the transfer of energy 
between different wave-numbers owing to the non-linear velocity term— 
become negligible. Thus the energy of Fourier components of very small 
wave-number is unchanged throughout the history of the homogeneous 
turbulence. In particular, this analysis shows that if the homogeneous 
turbulence is created with an anisotropic distribution of energy for the very 
small wave-number components of the turbulence, the anisotropy will 
persist throughout the decay. 

In general, anisotropy of the small wave-number components will escape 
notice in an experiment, because they normally contain a very small pro- 
portion of the total energy. However, there are some circumstances in 
which the permanence of the large-scale components emerges. For instance, 
the correlation between velocities at two widely separated points is strongly 
influenced by Fourier components of small wave-number and a lack of 
isotropy will be made apparent by a comparison of correlations over large 
space intervals taken in different directions. A lack of isotropy of the large 
eddies may also be expected to show up clearly in the final period of decay, 
because then all except the very largest eddies have decayed (1). It is the 
purpose of this note to present experimental evidence to show that a grid 
of bars in a uniform stream'‘creates homogeneous turbulence in which the 
components of small wave-number are anisotropic, and remain so through- 
out the decay as predicted theoretically. 


2. The spectrum at small wave-numbers 
If the general double-velocity correlation (in a field of homogeneous 
turbulence) between i- and j-components of the velocity at the points P 
—> 
and P’, where PP’ r, be written as Ri(r), u;u;, the spectrum tensor 
is defined as 


u 


rik) = — | | | Ri(r)e-**” dr(r), (2.1) 
where k is the wave-number vector, dr(r) is an element of volume at the 
position r, and the integration is over the whole of space. The correspond- 
ing Fourier transform relation is 
Ri(r) = | | Ti(k)e“*" dr(k). (2.2) 
The condition of incompressibility of the fluid requires I'(k) to satisfy the 
equations - ‘ i 
: k, Ti(k) = k; T%(k) = 0 (2.3) 


for all values of k. 


It can be shown that, as a consequence of the definition of T(k) and of 
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the continuity condition (2.3), the first non-zero term in the expansion 
of I'4(k) as a power series} in the components of k is 


kim ky U4 (2.4) 


imn* 
Here [i,,,, is a numerical fourth-order tensor which is required by the 
continuity conditions (2.3) to satisfy the relations 


ly ar eae Ti ant liane tT? 0. (2.5) 


im in} ijm imn mnt nim 
The dynamical equation for the rate of change of T%(k) with time then 
shows that I, is invariant throughout the decay. Hence, as stated 
already, the initial characteristics of the spectrum at small wave-numbers 
(e.g. anisotropy) will persist. The term (2.4) normally has very little 
influence on quantities like the total energy of the spectrum and the 
consequences of the invariance of T4,,,, only appear in rather special 
circumstances, two examples of which are discussed in the next two 


sections. 


3. The velocity correlation at large values of r 
We find from (2.1) and (2.4) that 

= ad 21... - =; | | rnt, Ri(r) dr(r). (3.1) 
This integral is invariant with time and will be isotropic only in the 
unlikely event that the large eddies are arranged isotropically at the 
instant of creation of the homogeneous field. The integrand is weighted 
strongly in favour of large values of r ( r|), so that any anisotropy of 
R(r) which exists may be expected to show up at large values of r, even 
during the initial period of decay when the turbulent energy is distributed 
with approximate spherical symmetry. 

From an experimental point of view, the best method of testing isotropy 
of Ri(r) is to choose i and j to correspond with velocity fluctuations 
parallel to the uniform stream in the wind-tunnel, and to vary the direction 
ofr. If x, is the downstream coordinate (with origin at the grid), a com- 
parison of R}(0,7r,0) and R}(0,0,7) (which should be equal in isotropic 
turbulence) provides a convenient test of symmetry about the 2,-axis. 
However, conventional grids consisting of two orthogonal sets of parallel 
bars are not likely to show very much departure from symmetry about 


] 


the x,-axis, even in the large eddies. A grid consisting of one set of parallel 


In reference (1) a pure imaginary term of the first degree in k,, was thought to be the 
t non-zero term of the series. However, since that paper was written J. Kampé de Fériet 
2) has derived an (exact) general form for the spectrum tensor, and consistency with this 


general form demands that the first-degree term in the series should vanish identically. 
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rods, 0-95 cm. in diameter with a spacing length 2-54 em., was therefore 
constructed and the above comparison was made. 

With a mean stream velocity of 620 cm./sec. the field was found to 
become statistically homogeneous at a distance 10 spacing lengths down- 
stream from the grid, at least in so far as the factors U,/Ug and E/E, both 
become unity, as shown in Fig. 1; U here denotes the mean stream velocity 
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E the turbulent energy, and the subscripts « and f refer to the conditions 
respectively behind and between rods of the grid. The energy of the 
turbulence is then distributed with approximate spherical symmetry, as 


9 


shown by the plot of u3/uz and u3/uz (where the x-axis is parallel to the 
rods). The isotropy does not, however, extend to the larger scales of 
motion, and the correlation coefficients 


Ri(0,7,0)/R1(0,0,0) and R}(0,0,7)/R1(0, 0, 0) 


are very appreciably different at large values of r. Fig. 2 shows the 


comparison of the two correlations at distances of 10, 20, 35, and 60 
spacing-lengths downstream. As is to be expected, the correlation between 
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e = velocities at two points with a relative displacement parallel to the rods is 
+ greater, at large values of r, than that for a relative displacement across 
o 4 the rods. 
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4. The final period of decay 
, mr ° ° “1s 1 ° . 
lhe second method used in this paper to exhibit the anisotropy of the 
> small wave-number components of the motion is to measure the directional 
the , distribution of turbulent energy in the final period of decay. When the 
60 


decay reaches an advanced stage, the velocity fluctuations become feeble 


and the motion is dominated by viscous forces. Under these circumstances 
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the variation of the spectrum with time is given by (1) pr 
I(k,¢) = T4(k, t))e-2**¢-0), (4.1) | 7 
When ¢—f, is very large, only the behaviour of I'i(k, t)) at small values of & sk 
is relevant. We have already seen that the expression (2.4) represents the ' 
spectrum at small values of & at all times, so that | af 
T(k, t) > k,, kh, Tenn C72 (4.2) | 
as t—t, > 00. In other words, the characteristics of the small wave-number 
components of the motion at the instant of creation of the homogeneous { 
field eventually determine almost the whole of the spectrum. 
The correlation tensor corresponding to the asymptotic spectrum tensor | 
(4.2) is found from (2.2) to be 
wd | )\4 Y . y 2 
¢i(r) - ao Cam — ™lexp| — _ ; (4.3) | 
OL 0 v(t ty) 8v(t— to) 
The mean squares of the three velocity components are obtained by 
putting r = 0; j 
ur uz uz (73/2)? 
i. I? 3 Siy(t—to] — 
mm 2mm 3mm 0 
The longitudinal correlation coefficient, usually denoted by f(r), is given by ? 
f(r) Ry(r, 0, 0) exp 4 | (4.5) 
; ur Sv(t—to) 
in view of the requirement of the continuity condition (2.5) that [},, = 0. | 
Thus the longitudinal correlation coefficient is independent of the direction | 
in which the correlation is taken, and has the form (4.5) whether the turbu- 
lence in the final period is isotropic or not. | 
Experimental confirmation of these predictions about turbulence in the — | 
final period of decay behind a conventional grid of two orthogonal rod 
systems has already been described (3). It is shown in that paper that 
measurements of u?, where the suffix 1 refers to the direction of the wind- + 
tunnel stream, vary as (¢—f,)~?, in agreement with (4.4), and that the, n 
longitudinal correlation coefficient, measured in the downstream direction, h 
has the shape given by (4.5). However, at the time when these experi- 
ments were made, the possibility of the turbulence not being isotropic in , ; 
the final period after being apparently isotropic earlier in the decay was 1 
not appreciated and measurements of w3 and uw} were not made. ' 
These measurements have since been undertaken with the same grid. t 
As was pointed out (3), the grid used, although appearing regular to the 
i 


eye, produced appreciable variations in uj in the plane normal to the 
stream at large distances downstream. These spatial variations were also 





Ek 
he 


LS 


or 


on 


yn, 
Ti- 
in 


yas 


she 


he 





ANISOTROPY OF THE SPECTRUM OF TURBULENCE 7 


present in uw and u3. However, the mean of a considerable number of 
measurements of each of u?, u3, and uz in the plane normal to the stream 
varied smoothly with distance downstream, and the averaged results are 
shown in Fig. 3. 

The ratios u?/uz and u?/u2 (which are equal within experimental error), 


after being close to unity in the initial period of decay, increase to about 
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1-5 far downstream. Also shown in the figure are calculations of dAj/dx,, 
made from measurements of Aj recorded in reference (3), where A, is defined 
by iF AR 9.9 (4.6) 
or r=0 
The asymptotic value of (U/v)dA}/dx, in the final period of decay (see (4.5)) 
is 4 and the curves show that this value is attained in the neighbourhood of 
the position at which u?/u2 approaches 1-5. The figure 1-5 (which may vary 
between 1-3 and 1-7 at particular positions across the flow) appears to be 
independent of the grid Reynolds number. Presumably, close behind the 
grid, the system of rods parallel to the 2,-axis produces large contributions 
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to uw? and w2 while the system parallel to the 2,-axis produces large 
contributions to uj and wu, the net result being that the contribution to w3 
from the largest eddies is permanently greater than the contributions to 
uz and w2. 
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By A. E. GREEN and W. ZERNA 
(King’s College, Newcastle-on-T yne) 
Received 28 June 1949] 


SUMMARY 


\ theory of thin elastic shells is developed as a first approximation without using 
the usual assumption that the normals to the middle surface of the undeformed shell 
remain normals in the deformed shell. General coordinates are used and the effect 
of shearing stresses is included so that the theory for plates as developed by Reissner 
is included as a special case. Vector and tensor notations are used throughout. 


1. Introduction 

In recent years a great number of papers have been published about the 

general theory of thin elastic shells. The methods used in these papers 

are varied and the theories are based on certain approximations which 

are often rather confusing and difficult to justify. Also there is some 

disagreement over the differences in the final forms of the results. For 

the usual theory of shells the following assumptions are made: 

1. The thickness ¢ of the shell is small compared with the other dimen- 

sions. Thus, if J is a characteristic length, for example, the smallest 
radius of curvature of the middle surface of the shell, thent+ 


A=t/L <1. (1.1) 


2. The displacements are sufficiently small for quantities of the second 
and higher orders to be neglected compared with those of the first 
order. 

3. The stress components normal to the middle surface are small com- 
pared with the other stress components and may be neglected in the 
stress-strain relations 

1. The normals of the undeformed middle surface become normals of 
the deformed middle surface. Sometimes the further assumption is 
made that the normals suffer no extension. 

The first assumption defines what is meant by the term ‘thin shell’ and 
may be regarded as satisfactory for many problems. The second assump- 
tion is always made in the classical theory of small displacements and 
small strains in elasticity and it ensures that the differential equations 


A is not to be confused with one of the elastic constants of Lamé, which is not used here. 


[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 
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remain linear. The third assumption may be reasonable, but the third 
and fourth assumptions taken together are much less satisfactory and 
are difficult to justify in a consistent way. Also it is difficult to estimate 
the effect of these assumptions on the orders of magnitude of certain 
terms which some writers retain and some reject. These points have been 
discussed in some papers, but no satisfactory conclusion has yet been 
reached. 

Using the above assumptions, general theories dealing with the equili- 
brium of thin shells have been developed by Love (1), Trefftz (2), 
Reissner (3), Reutter (4), Byrne (5), and Fadle (6) in which the lines of 
curvature of the middle surface are used as coordinates. Theories which 
are based on general coordinates are developed by Goldenweiser (7), 
Rabotnov (8), and Zerna (9). It is beyond the purpose of this short 
review to discuss all these papers in detail. They differ mainly in the 
methods used for the derivation of the fundamental equations and in the 
retention or neglect of certain terms which are connected with the initial 
assumptions. For example, Love (1), Reissner (3), and Goldenweiser (7) 
neglect certain terms, while Byrne (5) and Zerna (9) keep these terms and 
Byrne emphasizes that they should be kept. 

Papers which deserve special mention are those by Chien (10, 11, 12) 
and also a paper by Synge‘and Chien (13), in which, for the first time, a 
theory is developed which does not use the assumptions made by other 
writers. Chien’s theory seems to be remarkable for its complete generality. 
All quantities are developed in infinite series in terms of a coordinate x, 
through the thickness of the shell, and it is shown that it is theoretically 
possible to determine all the coefficients of the series. For practical 
calculations, however, it is necessary to approximate by considering only 
a finite number of terms of the series. It is here that some doubts must 
be expressed about certain of the approximations made by Chien. He 
appears to assume that the order of magnitudes of the various quantities 
can be assessed by expanding them as power series in A, and that differen- 
tiation with respect to coordinates in the surface of the shell does not 
alter the order of magnitude. Simple examples, however, show that these 
assumptions do not hold in many cases, for some of the quantities in 
question may depend on functions of the type exp(-_2/A"), where «x is a 
coordinate in the surface of the shell; these functions cannot be expanded 
in a power series in A and their orders of magnitude are altered when they 
are differentiated with respect to x. 

It does, in fact, appear to be difficult to decide how many terms in the 
series expansions in powers of x, should be retained. Also the order of the 





final differential equations which are obtained depends on the number of 
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terms retained in the power series, and this in turn affects the number 
of boundary conditions which it is necessary to impose on the edges of the 
shell. This introduces considerable complications if the order of the 
differential equations is too high and may make the theory too unwieldy 
for practical applications. 

\gain, Chien’s theory is developed from the ‘intrinsic’ point of view so 
that he avoids using the displacement vector. Although this is theoreti- 
cally satisfactory it has some disadvantages in practical applications 
where boundary conditions are often expressed in terms of displacements. 

[t appears that it is still desirable to attempt to develop a theory of 
thin shells which can be used for practical applications and which is, as 
far as possible, a systematic first approximation. At this point some 
guidance about the method to be used can be obtained from the theory 
of thin plates, since any theory of shells should contain the plate theory 
as a special case and since the plate theory is based on a reasonably 
consistent set of approximations. 

A theory of bending of thin plates which includes the effect of trans- 
verse shear deformations has been obtained by Reissner (14, 15, 16) by 
using Castigliano’s theorem of minimum energy, and the results are 
expressed in terms of certain weighted averages of the displacements and 
stresses. A similar method can be applied to the stretching of a plate. 
Green (17) has shown that some of the algebraic complications can be 
avoided, and the same results obtained, by using directly the stress-strain 
relations and the stress equations of equilibrium. It appears that some 
similar method to that used by Green for plates should be the most 
suitable to apply to shell theory, and such a method is developed in the 
present paper. There are, of course, additional difficulties in the theory 
of shells owing to the geometrical quantities which are involved. As well 
as using the methods of approximation which are similar to those used in 
plate theory, all terms of order A, or of higher order, are systematically 
neglected when compared with unity. Attention is confined to shells of 
constant thickness so that A is constant, but, apart from the additional 
algebra which is involved, there is no difficulty in extending the theory 
to shells of variable thickness provided that the derivatives of A with 
respect to coordinates in the middle surface of the shell are of the same 
order as A. 

The theory is developed with the help of vector and tensor calculus. 
Small Greek letters used as indices take the values 1, 2, 3 and small Latin 
indices take the values 1, 2. Unless otherwise stated the usual summation 
convention holds, so that summation over either of these ranges is signified 
by the repetition of an index. If the same index appears more than twice 
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in one term of an equation it is not summed. Partial differentiation is 
indicated by a comma and covariant differentiation of tensor quantities 
is indicated by a vertical line. Thus, for example, 


Oy 
a g,=> 
xX, F ’ 
a 
AyiB Ay—g—VXga,, 


€ 


Ye 2) 
a xB ly a “By I vy VB —I py me? 


where I3,, are the Christoffel symbols of the second kind and @, are 
independent parameters. The suffix in 6, does not indicate that these 
parameters are covariant since, in general, they have neither covariant 
nor contravariant properties. The differentials d6*, however, are contra- 
variant, and this is indicated by using the upper index. Also the Kronecker 
delta has the usual meaning 


i. = 1 (¢ =k, k not summed), 


2. Geometrical relations 

The points of the undeformed shell may be described by general curvi- 
linear coordinates x,. The surface which is determined by 2, = 0 is called 
the middle surface; 2, is the perpendicular distance of a point from the 
middle surface, and if ¢ is the thickness of the shell 

lt <a, < +t. (2.1) 

The faces of the shell are determined by the equations x, = +4t. The 
parameters x; represent any general coordinate system on the middle 
surface. 


Dimensionless coordinates may be introduced by putting 
~ 
0, : #,=—* (¢ not summed), 
: t 


where L, are constant characteristic lengths which are chosen in such a 


way that c<6.<d. 


u u i 
where the interval (c;,d;) is of unit length. Also, from (2.1) it is seen that 
6, lies between +}. 
The position vector of a point of the shell may be denoted by 
R* = r*¥+2,@,, 


where r* is the position vector of a point of the middle surface and e, is 
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the unit normal vector to the middle surface. Non-dimensional vectors 
R, r are now defined by the equations 


» 
R , = 
ee 


where L is a characteristic length which is the minimum radius of curva- 
ture of the middle surface, or the minimum of the lengths L;, whichever 
is the smaller. The points of the shell are now determined in non- 
dimensional form by R = r+aAd,e,. (2.2) 


[In the rest of the paper only R and r will be used, and these may be 
called the reduced position vectors. 

The expression (2.2) is the vectorial representation of a three-dimensional 
manifold with general coordinates 6@;,4,. The covariant base vectors are 
given by a, = R, = (6;—A6, b)e,, a, = R, = Ae;, (2.3) 
where e r, are the covariant base vectors of the middle surface, bf are 
the mixed components of the curvature tensor which is determined by 
the coefficients b;,. of the second fundamental form of the middle surface. 
Latin indices are raised or lowered with the help of g* and g;;, the contra- 
variant and covariant components respectively of the metric teasor of the 


middle surface. These are connected by the relations 


fix = ©", = 9g”, Jes = C.-C, = gg", 
Jia = C1 °C qq), Jn gik = dF, 
where Y Jik| = 911 J22—Iiz Gir: 


In addition 


e—gke, eo. = 9,.e’, e;, = T,e,+5;,e3, 
e.-e* — 5, e,-e, = 1, e'.. Ti,e’+bi.es, 
bi, = €g°@;, @; 3, in = Ce, = Ce; . (2.4) 


l rT’ 

A ae g ‘I rik? 
e. ble. eC, °C. | b 07 . | a eC, °C; x 
1 

>(Yri te T rk, Jik,r) 


The covariant components of the metric tensor of the manifold (2.2) are 





ML 5 y, a a Y ik 7 2A0, Diy T r°65 bi. bin, 
M ss a.-a 0) (2.5) 
Moe a,:* as A. 


If the determinant of the matrix m,, is denoted by m, then vm = ,/(g)m*, 
where m* = \(1—2A0, H-+-d262 K), (2.6) 


and where H = 3b” is the mean curvature, K = b}b3—b} 03 is the Gaussian 
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curvature. of the middle surface. The contravariant components of the 
metric tensor are found to be 


mm! = Mg Mgp, mm? = —mMy>. M33, mm** = m,, Mes, 
ms — Q, A*m** — 1, (2.7) 
and m‘* may be represented by infinite power series in AQ3, the first terms 
of which have the form 


mk — gik1 20, b+... (2.8) 
3. General state of stress 
The stress across the surface 6, constant is defined by stress vectors 


p, and these may be expressed in terms of the contravariant stress tensor 
s* by the relations 


i ae Ap. (3.1) 

Using (2.3) these become 
Po = — (o**e;,+-Ase5), (3.2) 
where ok — (8% —NO, bk)s™. (3.3) 
If Py, = V(mm™)p, (3.4) 


the equations of equilibrium. for an infinitesimal volume element cut out 
of the shell take the simple form 


D... = 0 (3.5) 
when volume forces are neglected. 

When discussing the equilibrium of thin shells it is more convenient to 
consider stress-resultants and stress-couples through the thickness of the 
shell, and the three-dimensional equations (3.5) are reduced to a two- 
dimensional form. Thus, firstly, equations (3.5) are integrated with 
respect to @, through the thickness of the shell and the resulting equations 
are expressed in terms of the components of the vectors with the help 


of the formulae (3.2) and the relation 


~ 1 ong 
, Soe 


NJ 06; , 


This process gives 


Nrk|,—AdE Q" AER’ + Pk = 0, ) 


z (3.6) 
N"b,,+-AQ"|, +AP3 = 0, J 
where the covariant differentiation means 
Ne, = NRL NLT Ne) oo 
Ol, = & +I J 
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The quantities V“*, Q' are given by the following integrals: 


Nik ( m*a'* dé, (3.8) 


t 


Qi = | m*s% dé, (3.9) 


and P*, R* are determined by the loads on the faces of the shell so that 


P° {he meee a R- 1m*@, 83%! 


(3.10) 


1° 


To obtain equations for the stress-couples, (3.5) is now multiplied by 6, 


and integrated through the thickness of the shell. This gives 


M* Qk —1ADEPr+ Re = 0, (3.11 a) 
Mb +AT"|, AGAR = 0, (3.11 b) 

where M*|, = MELTS M*+Ty, M*,) 
(3.12) 

Tr) = T 41%, 7, 
and M | m*0, 0% dO, (3.13) 
T’ ( m*0, s3* dé, (3.14) 
Q? — | m*s3 dé. (3.15) 


The Christoffel symbols in (3.7) and (3.12) are to be calculated from the 
metric of the middle surface, using the last formulae in (2.4). 

Equations (3.6) represent the conditions of equilibrium of the stress- 
resultants and equations (3.11) give the conditions of equilibrium of the 
stress-couples. The above derivation of these equations appears to be 
more straightforward than that which has been given by other writers. 

The quantities defined in (3.8) to (3.10) and (3.13) to (3.15) are of a 
purely mathematical character. It is not difficult, however, to interpret 
most of them physically. NN“ are related to physical stress-resultants 
and M* are related to the physical stress-couples; each forms a contra- 
variant tensor of the second order. Q‘ correspond to the shearing forces 
and form a contravariant vector. P* and R‘ are related to the components 
of the loads and the moments of the loads about axes in the middle surface. 


See e.g. W. Zerna (9). 
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4. General state of deformation and stress-strain relations 


The position vector of the deformed shell may be denoted by R* = LR 
and the displacement vector by V* = LV. The points of the deformed 
shell are now defined hy the non-dimensional vectors 


R = R+V. 
The covariant base vectors of the deformed shell are given by 
ay ay ry x? 


and the covariant components of the metric tensor of the deformed shell 


are M57. M34 a;-V,. T a, V; | 
Mis Ae;- V4 a; Vs, | (4.1) 
my) ) 7 
M33 = M33 2Ae,- V5. 


The covariant components of the strain tensor are now found to be 


Cag = 2(Mag—Mgg), (4.2) 


and using (4.1) these become 
€;, = 4(a;-V;.+a,-V,), 
€;3 = $(Ae,- V;+a;- V5), (4.3) 
€33 Ae; - V3. 
Assuming an isotropic homogeneous medium the stress-strain relations 
can be put in the form 


(1—2n)s** pl (1—2n)(m!*m** +-m*'m**) 4 2ym™*m"e,.4 en 
(1— 27)s*8 2um3| (l—7)m*%e,,-4 nme rl |, 
a** — 2um*me,,, 
(4.4) 


where is the shear modulus and 7 is Poisson’s ratio. 


5. Approximation 

Up to now all equations are exact in the sense of the classical theory 
of elasticity. No assumptions about the thickness of the shell and the state 
of deformation have been made (except, of course, that the deformations 
are small). The exact problem of the shell is governed by the equations 
of equilibrium (3.5), the stress-strain relations (4.4) which can be expressed 
in terms of the displacements with the help of (4.3), and by the boundary 
conditions on the faces and edges of the shell. In this form the problem, 
however, is for most practical cases too complicated and approximations 


must be introduced. 
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In developing the theory which is given here the following points have 
been kept in view: 

|. The theory is required to be a first approximation and regards stress- 
couples and stress-resultants as being of comparable importance. (This 
means that the approximation is beyond that of the so-called ‘membrane’ 
theory in which stress-couples are neglected.) 


» 


2. As far as possible the approximate theory should be developed from 
a consistent set of assumptions. The assumptions made here are 

a) the equations of equilibrium (3.5) will not be satisfied exactly but 
will be replaced by the equilibrium equations (3.6) and (3.11) for 
the stress-resultants and stress-couples. 
The thickness of the shell is assumed to be sufficiently small so that 
(1.1) holds, and all quantities of order A, or of higher order, are 
neglected compared with unity. 


With assumption (b) equations (4.3) become 


( s(e;-V,.+e,-V;,). 
“i3 3(Ae,- V; T e;- V5), (5.1) 
&33 Aes: V3. 


The terms containing A in (5.1) cannot be neglected since the relative order 
of magnitude of V; and V, is unknown. 
The displacement vector V may be expressed in terms of its covariant 
components bv the relation 
V ys e’ v,e*, 


and derivatives of V are given by 


V (v b Vz )e" +- (Ve 


dit 


! ra, \a3 
b’ v,)e%, 


| (5.2) 


vhere covariant differentiation is to be taken with respect to the middle 
surface. Thus 


, VW 
Uri Di Up 


while the component v3, which is perpendicular to the middle surface, is 
an invariant with respect to transformations of surface coordinates 6;, and 
therefore its covariant derivative is equal to its ordinary derivative, 1.e. 


v,;. Using (5.2) the strain components (5.1) become 


2b ik V3), 
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Again, using (2.7) and (2.8) and making the above assumption (6) about A, 
gain, g g ; 


the stress-strain relations (4.4) become 


9 
(1—2n)s HO 20a sHa)+ 20a a+ Fe oe | 
(5.4 


ag Qe(1— 5. 
A*(1 27)s*8 = sai os T 2ung”e,,; 


It is convenient to eliminate e.. from the first two sets of equations of 
33 | 

(5.4) in order to express the three stresses s‘* in terms of the three strains 

e;,,, and the stress s*°, Thus 


‘Ie vik o 7 Heo3! oie 
gtk pHttre, +A? 7 g*s*, (5.5) 


9 
4k : i 27 
W here Hikir gigk a gg 4 i] 


gi*g" (5.6) 
7 
is a tensor of the fourth order with the symmetric properties 


E ikir Rr ilr E ikrl Klrik 


The method to be followed in the rest of the paper is suggested by the 
corresponding method which has been used by Green (17) for plates. 

Expressions for the stress-resultants NV“ and stress-couples M“ are 
obtained from (5.5) by integrating both sides with respect to 6, after 
multiplication by m* and m*@, and using the approximation (b) about the 
order of A. Hence 


Ne—y_t gi*Q? = prAE*"(V,—b,, Us), 


i (5.7) 
Mi—— g*T? prAB*r(W, —b, We), 
where i= | 1d, = We= | 050,00, (5.8) 
U,= | r3d0, We= | 0505 d63, (5.9) 
+4 
z= | m*8, 83 dé, (5.10) 


The assumption (b) used in deriving equations (5.7) implies that 


+4 
v, 13, v0, d03, 5.11) 
x x3 ( 
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+4 
are of the same order in A and that | v, 63 d0, (n > 2) is not of greater 
rd 


order. Similarly 


[ vin G6, — | v4, a8, (5.12) 

= 
are of comparable order and | v; ;,63 d@ (nm > 2) is not of greater order. 
yf It should be noticed that no assumption is made about the relative orders 


1S 
of magnitude of | v; d6® and | v; ;, dé. 

:) With the above assumptions only displacement terms of the form (5.8) 
and (5.9) appear. If powers of A beyond the first are retained in the 
geometrical terms further displacement terms must be included. This 

9) would, however, take the theory beyond the state of a first approximation 
and is not required here.t+ 

It is now necessary to make some assumption about the distribution of 
the shear stresses s*’ through the thickness of the shell, and for this purpose 
| a power-series expansion in @, is appropriate. Since the present theory is 

1€ : ' , : a a : : . 
only concerned with a first approximation it is found, as in the theory of 
plates, that a quadratic expression in 6, must be assumed. If the boundary 

re ba Pe ; ‘ 
conditions (3.10) at the faces of the shell 0, ++ are then used the 

er 2 

expression for m*s* takes the form 

ie 

m*s3i — 3Qi(1—462)—1 Ri(1—1262)+ Pid,. (5.13) 
Using (5.13) the quantity 7" defined in (3.14) is 

7) o +a . 

1" ‘5 P*. (5.14) 
The third equation in (5.4) for the shear stresses s*‘ is now multiplied 
| by $m*(1—46%) and is integrated with respect to @,. This process avoids 

3) ee 3 : ; . 
the introduction of unknown values of the displacement at the surfaces 
of the shell and, with the assumption about A, gives 

».9 f ’) , ;, y . y 4 = 

) A($Q'—1R*) pg" (AV, + 12M). (5.15) 

10) where A 3 (1—462)v, d9, (5.16) 

the remaining quantities being defined already in (3.9), (3.10), and (5.8). 
Byrne (5) retains quantities up to order A’ in the geometrical terms but, as Reissner (18) 
11) ind other writers have already mentioned, this does not appear to have much meaning since 


some terms of order A have already been neglected in other parts of his theory. 
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From the second equation in (5.4) it follows in the same way that 


1 21) og A 1 grt | (1—463)v,,, de (5.17) 


= 
24u(1 ”) 8] 1 


where S3 — 3 | m*(1—46§)s* das. (5.18) 


Again, multiplying the second equation in (5.4) by m*6,(1—463) and 


integrating gives, after some simplification, 


” id . oe -2n =: - 
g”™ | 0,(1—46§)v,., de 4 S3 (5.19) 
Zl—yn J : 4 p(l—n) 


where S3 [ m*0,(1—463)s*3 dé. (5.20) 


Before reducing the expressions (5.7) for the stress-resultants and stress- 
couples to their final forms it is necessary to obtain expressions for Q? and 
T°’. Using equations (3.11 b), and (5.14), Q* becomes 

AQ? = Mb, +-BAP"|,+AR3. (5.21) 

To obtain a corresponding expression for 7% the equation of equilibrium 

(3.5) is multiplied by 6% e, and is integrated with respect to 6,. With (5.13) 


and the second of equations (3.6) this gives 


ATS LN) +41", + LAP3+ AAR’ |, (5.22) 


where Lik | m*03 ot dO, 


[f now (5.17) and (5.19) are substituted in (5.7), and if the assumption 
about A is used, the following equations are obtained for N*, M*: 





N%& = prAH*"*(V,—b,, Vy) +X", (5.23) 
M* = prAR“Wp, +Y*, (5.24) 
he =) = 
where Xtk “” gk be oi b,, E83 |, 
] 7) 4n 
Vik A*n gikps AC “2) 5 Rikirgs ; 
I—y|' 24n =” . 


Because of the assumptions which have been made, the first groups of 
terms on the right-hand sides of (5.23) and (5.24) are of the same order 
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of magnitude and the quantities X“ and Y“ can, at most, be of this same 
order, and may be of smaller order. The stress-resultants N“ and stress- 
couples M‘ are therefore of comparable order. Also, the quantities S* 
and 7° are, in general, of comparable orders of magnitude, as are also 
S? and Q%. It can then | 


e shown that, whatever the relative orders of 
magnitude of 7° and Q%, the terms in S? and S? in X“* and Y“ respec- 
tively may be omitted. Then if AQ*, AT? are substituted from (5.21) and 
(5.22) in (5.23) and (5.24) respectively, and, if it is remembered that in 
general, AL** can be neglected compared with M“, the stress-resultants 
and stress-couples can now be expressed in terms of the five unknown 
quantities V,, W;. Thus 


l 


A? , - ; . i es 
, l 4 h° 13" r) pAktkr(V, bV3), (5.25) 
| 
° } 3 Ll Pr Yikir |Z 5 OF 
M* — Tie g**(P3+2R"|,) pA Wy. (5.26) 
Ui 
The shearing forces Q* have already been expressed in terms of the 


displacements V; and W; in (5.15) and the expressions for these forces are 


repeated here fo completen SS 


6Qi_1R 5 7(ANa + 12M). (5.27) 


4 


The shell problem is now governed by the five differential equations (3.6) 
and (3.11 a) and the equations (5.25) to (5.27). 

In the speci il case of the flat plate the above theory reduces to that given 
by Reissner (14, 15, 16) as far as bending is concerned. It should be 
noticed that Reissner’s theory for the bending of a plate requires that 
three boundary conditions are to be satisfied at the edges of a plate in 
contrast to two boundary conditions in the classical bending theory. In 
dition, the stretching theory for plates requires two boundary conditions, 
making five in all. Similarly, five boundary conditions are required in the 
present theory of shells, in contrast to four in the usual theory. 

When the shearing forces in (5.27) are neglected, and orthogonal 


curvilinear coordinates are used, a theory is obtained which does not 


contain certain terms which are introduced by some other writers, but, 
is mentioned earlier, the inclusion of these extra terms only appears to 
be necessary when proceeding beyond a first approximation. The justifica- 
ion for neglecting the shearing forces is not obvious but there may be 
many examples where this is satisfactory. The inclusion of shearing forces 


is likely to be of importance mainly for ‘edge’ effects. 
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ON THE PLANE STRESS-DISTRIBUTION IN AN 
INFINITE PLATE WITH A RIM-STIFFENED 
ELLIPTICAL OPENING 
By A. A. WELLS (Department of Engineering, University of Cambridge) 
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SUMMARY 


\ solution is given for the decrease of stress concentration factor at an elliptical 
opening in a plate under uniform plane stress at infinity when a compact stiffening 
rim is applied to the opening. We begin with the known solution for a plain elliptical 
opening, on to which is superimposed a stress distribution satisfying the boundary 
conditions between rim and plate and vanishing at infinity, and the distributions 
of cross-section of the rim corresponding to these are derived. Particular solutions 
are obtained for hydrostatic tension and pure shear in two directions at infinity, and 
from these all other cases may be solved. For ratios of major to minor axis which 
are not too great, and for rims of moderate cross-section, the cross-sections corre- 
sponding to the given solutions vary little round the opening. 


Symbols 
x, 8, elliptical coordinates. 
a, b, semi-major and semi-minor axes of elliptical opening. 

c, semi-interfocal distance for elliptical coordinate system. 

B, «8, stresses in elliptical coordinates. 
n, any positive or negative integer. 
P,, Q,, arbitrary constants. 

A, area of rim cross-section at any point. 

r, radius of opening at any point. 

t, plate thickness. 

y, A/rt. 

v, Poisson’s ratio. 
K,Y,F;,C, constants. 
Previous work 

The reduction of stress concentration due to a rim homogeneously 

attached to the periphery of a circular opening has received some attention 
in the literature (1, 2, 3). The rim is assumed to be of uniform cross- 
sectional area A, placed compactly at the edge of the opening so that 
the stress in it is uniformly distributed over the cross-section. If the 
radius of the opening is r, and the plate thickness is ¢, then the parameter 
y — A/rt defines the influence of the rim on the stress concentration factor. 
If the radial and circumferential stresses in the plate, at a point where 
the polar coordinates are r, 0, are 77, 66 and the shear stress is 7, from 


[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 





24 A. A. WELLS 


considerations of strain equality and equilibrium the boundary conditions 
at the periphery of the opening become 
06 = 7rr(1+-vy)/y, r0 drr/dé, (1) 
where v is Poisson’s ratio. At infinity all combinations of principal stresses 
may be obtained by superposition of mult iples of the solutions embodying 
hydrostatic tension, for which 
1? = 0 = I, rd = 0, (2) 
and pure shear, for which 
66 cos 26, rr —cos 26, rd —sin 26. (3) 
The stresses in the plate at the periphery of the opening, and the rim 
stress 4@, corresponding to these cases are, for hydrostatic tension at 
infinity, 7; , = , . 
: 60 2(1-4 vy) F,, rr 2yt (9 1) 
00, = 2F,, ro) = 0, (4) 
where F;, is defined as 1/{1+-y(j-+-v)}, and for pure shear at infinity 
y 2A DP , iy - ne D i" , ¢ 
66 4cos 26(1- vy)F;, r7 4 cos 20 yF; (7 3) 
66, 4 cos 26 F,, r0 8 sin 26 yF;. (5) 
The circumferential plate stress 60 at the periphery of the opening is 
always greater than the rim stress 06,, and the stress concentration factor 
is determined either by 6@ or r@, which assumes importance when y is large. 
For hydrostatic tension there is no limit to the economical value of y, but 
for pure shear the optimum is reached when y = 1/(4—v) on the maximum 
shear criterion, and the stress concentration factor is reduced from 4 to 16/7. 
All other cases fall between these two extremes. 


Boundary conditions: elliptical opening 

To maintain continuity of nomenclature the rim parameter for the 
elliptical opening is defined as A/rt, where 7 is now the radius of curvature 
of the opening at each point under consideration (Fig. 1). In the general 
case neither A nor r will now be constant round the periphery. As before, 


there is the condition, using the notation adopted by Inglis (4, 5), 


BB va(1-+-vy)/y, 
' xO 
which can be restated as y = - (6) 
BB—vacx 
For shear equilibrium along the periphery the equation now becomes 
dp di Aw 7 n 
: (ca) vB, (7) 


ds dB : 








wher 


coor 


By ¢ 
disp! 


gene 


elli 
de} 
of : 


in 


ag 











L) 





coordinates this is rewritten 


me ) 
6aa sin 26 


general expressions for dist1 


PLANE STRESS-DISTRIBUTION IN AN 








INFINITE PLATE 


25 


where s is the distance measured along the periphery. In elliptical 


vB sinh 2x». (8) 


By expressing the equations of equilibrium for plane stress in terms of 


displacements, and introducing the complex variable, Inglis has obtained 


Q 


ibutions of stresses aa, 68, and af in an 
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(n Q) 
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p> -p|>— S¥ 
o~ = 














3)P+-e°* sin(n- 1)B] 


p€——_——— [\) 
| a 


1 terms of infinite series, the terms of which 
nteger n and arbitrary constants P, and Q,,. Substitution 


nto equation (8) yields the following equation 


3)B+-e-*sin(n 5)B| 
1)B-++-e* sin(n—3)B | 


1)8+-e* sin(n— 1)p | 


3X sin(n 1)B+-e—* sin(n 3)p 


1)B+-e*sin(n—1)B] 0 


3)B8+-e3* sin(n- 1)8 


I 


(9) 


of expansion into a series, and the same is true whether the 
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functions of B in ac, xB in the general expressions are cos, sin or sin, —cos 
respectively. 

At infinity three sets of boundary conditions must be considered, to 
cover all principal stress combinations. These are 


1. Hydrostatic tension: 


vax = BB = 1, xB = 0. 
2. Pure shear: 


xx = —BB =cos2B, of = —sin 26. 


a BB =sin2B, of = —cos2P. (10) 


Subsidiary stress expressions 

Equation (8) is automatically satisfied by the stress expressions repre- 
senting no rim. Calling these stresses at the periphery aa, BB, “B, 
(<x, = 88, = 0), and those of any general expression, vanishing at infinity 
and satisfying equation (8), ads, BBo. xBo, by the principle of superposition 


~s 


va = Kit, BR=fB,+KBB, of = Kap, (11) 


where AK is a constant, represents a further solution of equation (8). 
Equation (6) for y now becomes 


y = K&x,/(BB,-+ KBB.—vK xa), (12) 


from which it is evident that 


TH ~B SB a || ry 2 ") 


y 1+-vy q XX: 
~ a Br 
while for the rim BB, a, /|1 uf - PB: 4 ‘| 
XX 


Subsidiary stress expressions for superimposition on those representing 
no rim, vanishing at infinity and satisfying equation (9), must therefore 
be found, with the additional qualification that a rational value for y in 
equation (12) be the result. 

By elimination of values of » giving finite stresses at infinity and 
unsuitable distributions of stress at « v9 in the general stress expressions 
the number of terms in the series is considerably narrowed. Suitable 
expressions are obtained by putting Q_, = 1 for hydrostatic tension and 
P, 1, 2, —2e*% for both cases of pure shear. Substituting these 
constants into the general stress expressions the following values of the 


stresses at a vy are obtained: 
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ay BB. 202sinh2a), af. = —2C*sin 28. (14 
2. Pure shear: 
vx, cosech 2 Xo 8 35 | (2 cosh Zao + sinh 2a) 2 cos 2B | 
4C?(e-2% — cos 28), 
vB. = 40? sin 28(e»—cos 28). (15) 


3. Pure shear at Ltr: 


a ¢ »Q ¢ 
vx, cosech 2x, BB,/| (2 cosh 


2 4 ('2(1 —e2 


XYVo 


where C (cosh 2~,— cos 28). 
Results 
BB, X 


The solutions for ax,, ; 
for a 


l. Hydrostatic te NSION . 


—-) 


2 cos 28] = 4C*sin 28, 


(16) 


2a,+sinh 2a) 


‘© cos 28)(1—e-?* cos 28), 


are well known, so that the final solutions 


v, from equations (13) become 


xO 38 a ‘ 3 vB a Doin 9 
2F. C sinh 2ap, 2F, C'sin 28, 
7] 1+ vy Y 
bp 2F.Csinh2a, (j 1). (17) 
> Pure shear: 
2D 
XX é 7 y ‘ y 
r 2F. C(1—e? cos 28), 
Yy l 7- VY , 
vB oF pe sin 2B(e*° cos 28) 
y sinh 2x, ; 
3B 2F.C(1—e?* cos 28) (j m). (18) 
3. Pure shear at bar: 
¥08 of) 
mnt 2F, Ce? sin 2B, 
y l+-vy , 
Me} oF.C! 1 — 2% eos 28)(1—e-?% cos 28) 
y - sinh 2a, 
3B 2F.Ce?sin2B (9 m), (19) 
2 cosh 2 sinh 2a 2 cos 2B 
where m “0 = : (20) 


These expressions reduce to 
becomes infinite. If y at fb 


(12) yields the following relation 


Y Yo 


sinh 2a, 
those for the circular opening when ap 


z be termed y,, substituting in equation 


>. aa 
XX PP 4-9 
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where a%5.9, BP1.9, BBs.9 denote their values at « = a, and 8 = jz. This 
reduces to . 
cosh 2 Xo 
Y Yo ) 
cosh 2a4,—cos 26} 1-+-y,(1-+v)} 
for the three cases dealt with. Superposition is therefore straightforward. 
Again the relation is y = Y_ for a» equal to infinity. 
To determine the stresses in the interior of the plate, K in equation (11) is 
given by 3 
K PP 1-0 Yo : (23) 


- ye) 
x%9.9 1+-Yo(—PP2.9/%%2.9+V) 


i 


At a ¥p Maximum values of «8 for cases 1 and 2, and ff for case 3, occur 
away from the major and minor axes of the opening. By substituting y, 
Q 


for y and differentiating with respect to 8, the values of tanf at which 


the maxima occur are given by 


L. Hydrostatic tension (xB): 

- (cosh? 2a ae 

tan?8 — % 0 
cosh 2a,)+ 1 

2. Pure shear (af): 

sinh 2«,[./(9-+- 2e2% sinh 2v9)—3] 


tan*8 ~ 
(cosh 2a,)-+ 1)(e2%+-1) 


3. Pure shear at hor (8 


Js 


[Se) 


cosh 2a] ./{9(1— Y)?+- tanh? 209(¥? cosh? 2a)—1)}+-3(] Y)| 


tan*f —— 
(cosh 2a,+ 1)(Y cosh 2a,+ 1) 
; l+v 
where ) Jo__ (24) 
1+ YoU v) 
When «a, becomes infinite these expressions give tanf 1, which 


agrees with the result for the circular opening. The remaining stresses at 
the periphery have their maxima at 8 = 0 or $7. It is important to note 
that the stresses to which equations (24) refer may be those determining 
the stress concentration factor at the opening if the amount of rim 
stiffening is large. 

Equation (22) gives different distributions of y on the periphery of the 
opening for the various values, both of b/a and yo. The area of cross- 
section, A, given by yrt, may be expressed as a fraction of Ag, that at jz, 
and this fraction has been plotted against 8 for b/a = 2/3 and 1/2 with 
various values of y, in Fig. 2. When yp, is small A tends to be larger at 
the minor axis than at the major, but this tendency is reversed as Yo 


increases. 
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Application of results 


The solutions assume distributions of rim cross-section at the periphery 


which vary slightly and may not therefore be wholly practicable. It would 


seem reasonable, however, to invoke St. Venant’s principle, by applying 


ie 


solutions to isolated pon 
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dictated by the design at these points. For this purpose equations (17), 


Yn 


LS), and (19) would be used as they stand. 


Given the values A, a, b, and t, y may be calculated as follows: 


If A 


h ad 


At the major axis 


9a ») 
Av2sinh 2a, = 
y . (25) 
(cosh 2a)— cos 2B)! 


stress values at the periphery at the major and minor axes 


ay be simplified and become 








30 A. A. WELLS 


At the minor axis 


l A 
m 2+h), y = k—. 26 
i Y at St 
1. For hydrostatic tension. 
At the major axis 
= 2(1+-vy) 


BB 


At the minor axis 


) 
bo 
=~ 
— 
= 
= 
— 


Ro 


ry 


; I+y(1-+v) 


2. For pure shear. 


At the major axis 


= l l+-vy 
3B 2;1+ : R 
PI ( ar +-y(m--v) 


At the minor axis 
- 1k tas 
6 2(1-+-k)(1 vy) (28) 
1-+-y(m-+-v) 
3. For pure shear at har. 
At the major axis 


B 2y(1+k)? 


k{l+-y(2k+1+)}° 
At the minor axis 

“ 2y(1-+-k)* 

Xf 


k+y{2+k(1+-y)}? 


Ws 


(29) 


and for the distribution of y round the opening specified by the solutions 


At the major axis 


A 2 \3 9 
. ( 2) 2 ; (30) 
A, \! -B) 2+ (1—k*)yo(1+r) 


where Ay 


Conclusions 


For an elliptical opening stiffened by means of a homogeneous rim an 


exact and rational solution is obtained when the area of cross-section of 


the rim follows a distribution round the opening which depends both on 


the eccentricity of the ellipse and the amount of stiffening. For values of 


these variables to be expected in practice the distribution does not depart 
far from uniformity. In practice also, calculations may be carried out for 
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any particular point by use of the value of rim parameter there and 
invoking St. Venant’s principle. By this means it is probable that the 
solutions presented may be applied, with only a nominal loss of accuracy, 
to elliptical openings having a uniform cross-section of rim. 

Except for the case of hydrostatic tension there is a maximum limit to 
the economical cross-sectional area of stiffening which may be applied, 
since reduction of the stress in the plate tangential to the opening is only 
achieved at the expense of boundary shear stress in the plate, necessary 
to build up load in the rim. This limit depends on the state of stress in 
the plate at infinity, and is determined qualitatively by equating maximum 
tangential and twice the shear stresses at the periphery of the opening. 

For the particular case of the circular opening all the solutions presented 
revert to those obtained by Timoshenko, Beskin, and others. 
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SUMMARY Ky 
By introducing the concept of equivalent elastic systems, it is shown how the L 
forces in the members of multiply-redundant pin-jointed plane frameworks may be | L’ 
obtained directly as the result of simple arithmetical operations only. The method P 
of analysis proposed is divided into two distinct parts of which one is an expression y 
of the elastic properties of the framework, valid for all conditions of loading, while 
the other takes into account the actual loading conditions in any given case. A 
A numerical application to the case of a rectangular-panel truss is given. 
: w 
1. Introduction ; 
ft 
[HE analysis of pin-jointed redundant plane frameworks by any of the ; 
classical methods involves the preparation of as many elastic equations 
as there are redundancies, and the simultaneous solution of these equations. 34 


In all but the simplest cases this means that a calculating machine is | 
needed to cope with the number of significant figures required if important : 
errors are to be avoided. 


In 
In the method described in the present paper there are no simultaneous loadi 
equations to be solved; only simple arithmetical operations are used, and 
an accurate solution is obtainable by the use of a slide rule only. The 
method is particularly suitable for use when several loading conditions 
have to be investigated. 
2. Notation 
The following notation will be used in this paper: fi 
A denotes cross-sectional area of a bar. 
a denotes cross-sectional area of an additional bar forming part of an 
equivalent bar. 
B (Ly3)° | (Lg4)? , (L44)°Lo3_ 
EA,, HAy EA,, 
B, (Lis4)” ; (Log) (Lo3)"Ly4_ 
EA, EAw EA, 
C denotes the elastic constant for a panel or group of panels. 
D_ denotes the relative axial displacement of cut ends in bar 23 pro a te 
duced by a unit tension applied in the line of the cut bar 23. a 


{Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 





















PIN-JOINTED 





REDUNDANT FRAMEWORKS 33 


D, denotes the relative axial displacement of cut ends in bar 14 pro- 
duced by a unit tension applied in the line of the cut bar 14. 

E denotes Young’s modulus. 

F denotes the force in a bar. 


| denotes the force in an additional bar. 


q kK B/¢ 
K, B,/¢ 
Kk 1—K 
K, l—K, 


L denotes the length of bar; force in a bar. 
ee | L’ denotes the extensibility (1/A E£) of a bar. 
id P denotes the external force applied to a framework. 
S’ denotes the extensibility of an equivalent bar. 
A denotes the relative axial displacement of two ends of a bar due to 
force F in the bar. 


w Li4/D, Log or D53/DL 44. 

















f f as suffix, refers to equivalent bar. 

ie . . e,° . ° . 

1, 2, as suffixes, refer to the two additional bars forming part of an 
lS ° 

equivalent bar. 

34, 23, etc., as suffixes, refer to bars 34, 23, ete. 
IS 
it 3. The single redundant panel 

In the general case of any redundant panel not subject to external 
IS loading, such as 1234 (Fig. 1), the internal forces in the bars produced by 
e 
oi 2 ~L 44 Y 
a 

; 4 
(b) 
‘1G. 1. Relative internal forces in bars of a redundant panel not subject to 


xternal loading. 


. force in any one member are related to the lengths of the bars as indi- 
cated in Table 1 (a), in which the positive sign is used for tension. This 


9 


D 
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can be shown to be so by considering bar 24 to be removed and replaced 
by a force + L,,; this will create forces in the other bars having the values 
shown in the table. It will be noted that the force in bar 23 is — L,,, and 














that in bar 14 is —Z,,; otherwise forces are proportional to the lengths 
of the bars. (These are the bars which, in a multi-panel structure, are 
common to adjacent panels, i.e. the vertical bars in the case of trusses or 
the horizontal bars in the case of towers.) In the special case of a rectangu- 
lar panel, all forces are proportional directly to the lengths of the bars. 
TABLE | 
Bar | 12 | 34S 13 ie. | -*+ | om 
(a) En ‘ 7m bie i - 
(by KL, | +K’Ly K’Lis + KLng KL phe Kin.: 
(c) £% | 40h. | +hike he | he~Ske K, Los 
(dy in | —nile, | +h | +0 has oP Li, (liad 1 
(a) Relative internal forces in bars of panel (no external loading). 
(6) Internal forces in bars of panel (external loading, case I). 
(c) Internal forces in bars of panel (external loading, case II). 
(d) Induced forces in bars of panel. FI 


If the bar 14 be cut, and equal and opposite forces + L,, be applied in 
the direction of the cut bar, it may be shown that the relative axial 
displacement of the cut ends is 


C 
Pate ] 
} ay ( ) 
: : : ae ae | a fs 
in which ( wat (Lo)? EA., +(L,,)* BA,’ 


the summation referring to chord members and diagonals. The quantity C 

will be called the elastic constant for the panel. On the other hand, if a unit 

force is applied, instead of a force L,3, the relative axial displacement of 
the cut ends is 0 

. » 

oe ” 

The effect of external loading will now be considered. First, if the 

panel 1234, supported as shown in Fig. 2, is acted upon by a horizontal 

force P of amount L,, at the joint 1, the other external forces required 

for equilibrium are as shown in the figure. It may then be shown that if 

p= (L43)° (Liq)? | (L444)? Los 


: BAys | BAgy © BAY? { 


[ = ae and K’=1-—K, 



















{elative internal forces in bars of a redundant panel due to external 
loading on the left (case I). 
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Fic. 3. Relative internal forces in bars of a redundant panel due to external 
loading on the right (case IT). 
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the internal forces in the bars due to the external loading are as given in 
Table 1(b). In the case of rectangular panels, a simplification is possible 
since L,, = L,3 so that the force in bar 14 may be written A’ L,, instead 
of (L,,—KL,3). 

Secondly, if the same panel 1234, supported as shown in Fig. 3, is acted 











v 
Y 
r 
+ A + J 
Fic. 4. 


upon by a horizontal force P of amount L,, at joint 2, the internal forces 
in the bars are as given in Table 1 (c). Here 


K, a, kK; 1—K, 
(Lg,)® , (Leg)® , (Leg)*L 

and B a a 4 a, 
' BA, ‘ BAy EF Ags 


For rectangular panels, the force in bar 23 may be written A; L,, instead 
of (L,,;—K,L4,). 

All bars in a panel may be replaced by one equivalent bar having the 
same elastic properties as the bars in the panel, considered as a group. 
Thus, the single panel illustrated in Fig. 1a may be considered as com- 
posed of the two parts shown in Fig. 4. If equal and opposite unit forces 
are applied at the joints 1 and 4 (Fig. 46) in the direction of the missing 
bar the relative axial displacement of the joints is 


C l 
-— 4 = D—Ly. (3) 
(Lg3)? EA, ; se 
It follows that to produce a unit relative axial displacement the force 
required is | 


=~ (4) 
D, Liy 
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To produce a unit relative axial displacement of the ends of the bar 14 


itself (Fig. 4a), the force required is 


l ” 
ae (5) 

| 
so that, for the complete panel (Fig. 1a), the force needed to produce a 


unit relative axial displacement of the joints 1 and 4 is 


. D, 


. a . — (6) 
Lin Lig Li4(D,—Ly,4) 


Elastically, a single bar between joints 1 and 4 is equivalent to the six bars 
in the panel pro\ ided that 

| D, 

S14 Li (D, Ly4) 


where S;, is the L/EA value for the equivalent bar. Hence, 


D,—L, “ 
S14 Lisl D, “), (7) 


If this equivalent bar is assumed to have the same length as the original 


bar 14, its cross-sectional area A, must be given by 


D, 
A, Auly, val (8) 


Since the term in brackets is greater than unity, the equivalent bar may 
always be considered as made up of the original bar of area A,, together 


with another bar of the same length and of area a such that 


’ L. 
? A : 14 . 9g 
‘= Aulp, ) si 


the two bars acting in parallel between joints 1 and 4. 

The effect of a force acting in the direction of such an equivalent bar 
may then be regarded as divided between the original bar and the addi- 
tional bar in the ratio of their areas (Fig. 5). The force in the additional 
bar induces a set of internal forces in the other bars in the panel in accor- 
dance with Table 1 (a) and these will be called the induced forces. Thus, 
suppose a pull P to be applied to the panel at joints 1 and 4 in the direction 
of bar 14 (Fig. 6). Then the forces in the bars of the panel may be obtained 


as (a) a force in bar 14 of amount 


P(D,—L44) 


> 10 
7 (10) 
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and (b) a set of forces in all other bars induced by a force in the direction 14 
of amount a 
I Ly4 


1] 





The force 
































' 
sponds to the case in which bar 14 is cut and has its cut sections forced 
(2) 
-+— I 
a =— 
‘ P , Ds ‘ 
4 (p= 
b) 
(0) 
f, bor : 
=< - tod ' 
Ae j 
Vd. 
—_—_—— ——; gs | 4 
7 \ 
4 (5 ) 
Fic. 5. Equivalent bar (a) considered as original bar plus additional bar (6). 
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Fic. 6. Pull on bar 14 tends to separate joints 1 and 4, 


apart, i.e. to the case of a compressive force PL,,/D, in the bar. By 
proportion from Table 1 (a) the forces induced in all the other bars may 
be written down as in Table 1 (d), in which w + Ly4/D, Los. 
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These forces are induced from left to right across the panel, but if P 
acts on joints 2 and 3 in the direction of the bar 23, the forces are induced 
from right to left. The values in Table 1(d) still apply to this case 
provided that w L5,/DL,,, where D = C/L#,. In this table the force 
juoted in brackets is the force in the additional bar with the sign reversed 


for the reason given. 


4, Multi-panel trusses 

In the case of a truss composed of a series of redundant panels, the end 
panel on the left may be reduced in the manner described above to an 
equivalent bar for which the S;, value may be readily determined. If this 
value is then used in place of the L,, value in the second panel, a second 
S3, value may be found which will be the equivalent bar replacing the 
first two panels. In this way an equivalent bar (of value S3,) may be 
found which can be used in replacement of any number of panels to the 
left. Similarly, another equivalent bar (of value S;,) may be found which 
can be used in replacement of any number of panels to the right. The 
elastic constant for any such group of panels may be found as for a single 
panel. In doing this, a panel containing an equivalent bar is considered 
and the S’ value for the equivalent bar is used in place of the L’ value for 
the actual member there. 

Any multi-panel truss may therefore, in general, be replaced by the 
chord members and diagonals of any one panel, together with two equiva- 
lent bars (of values Sj, and S;,) replacing all panels to the left and right 
respectively. Such a system is called an equivalent panel. It follows that 
as many equivalent panels may be considered in any particular case as 
there are panels in the truss. The elastic constants of such equivalent 
panels are found in a manner similar to that which was described above 
for the single panel 

The forces acting at the joints of any equivalent panel are as follows: 
(a) a vertical force at joint 3/joint 4 equal to the algebraic sum of all the 
vertical forces to the right/left of the equivalent panel, and (6) two equal 
and opposite horizontal forces at joints 2 and 3/joints 1 and 4 exercising 
the same moment as that of all the forces to the right/left of the equivalent 
panel. These forces are provided automatically by the interaction between 
the panel considered and others on either side of it; such forces will be 
called the conjugate forces for the panel. They correspond to the external 
loadings considered in the case of the single panel. Hence, values of B, B,, 
K, K,, K', and K; may be found for any given equivalent panel. 

The conjugate forces acting on an equivalent panel produce forces in 
the bars in accordance with Tables 1 (b) and 1 (c), and these are called the 
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primary forces. They differ from the actual total forces by amounts which 
represent forces induced from other panels. 

In any equivalent bar in an equivalent panel, a part of the primary 
force induces a set of internal forces in the bars of the next panel. Thus, 
if P is the primary force in equivalent bar 23 of the first panel on the left, 
the induced force in bar 14 of the second panel is, from (11), P(—L,,/D,), 
where L;, relates to bar 14 in the second panel and JD, relates to the group 
of panels to the right (including the second panel). The induced forces in 









































I I mm Vv 
a e 
40! 
Ln, 
A Nl + 
v Vv v 
g! 16" 247 
ae 4 @ 30° = 120’ > 
1 2 
Ne fo os for hor af 
ere right / 
4 A 3 


Fic. 7. Four-panel truss with four redundant bars. 


the other bars of the panel are related as in Table 1(d). Usually, the 
value w referred to in this table is very small, which means that the 
induced forces are relatively very small and therefore effectively vanish 
after being induced from any one panel to the next. In this way forces 
are induced from left to right from panel to panel. 

Similarly, if P is the primary force in the equivalent bar 14 of the last 
panel on the right, the induced force in bar 23 of the last-but-one panel 
is P(—L,,/D), where Ly, relates to the last-but-one panel and D relates 
to the group of panels to the left (including the last-but-one panel). The 
induced forces in the other bars are related as in Table 1 (d). Thus, forces 
are induced from right to left from panel to panel. In the case of a 
symmetrical frame loaded symmetrically, these two sets of induced forces 


are similar and do not require to be separately evaluated. 
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The total force in any bar is the sum of the primary force and forces 
induced from other panels. 
5. Method of analysis 

This is divided into-two distinct stages which will be called the structural 

and the stress analyses, respectively. 

The structural analysis involves the following two steps: 

1. Evaluation of the equivalent bars for successive groupings of panels, 
working from left to right and also from right to left if the truss is 
unsymmetrical. In the case of a symmetrical truss the two sets of 
values will be similar. 

2. Evaluation of the constants C, B, B,, K, K’, K,, and K{ for each of 

the equivalent panels. 

The structural analysis, which is an expression of the elastic properties 

of the truss, is valid for all conditions of loading. 

In the stress analysis the forces in the bars due to any particular loading 

condition are determined according to the following steps: 

3. Evaluation of the conjugate forces acting on each panel. 

4. Determination of the primary forces in all bars from Tables 1 (b) 
and 1 (c). 

5. Determination of all the induced forces working from left to right 
and also from right to left. 

6. Summation of the primary forces and the induced forces to obtain 
the total forces in the bars. 

6. Example 

To illustrate the method of analysis using equivalent elastic systems, 
the case of a four-panel truss having four redundant bars (Fig. 7) will be 
considered. This example has been analysed elsewhere (1) by the two 
classical methods of virtual work and least work, and also by the Williot 
strain method, and the solutions reduced in each case as far as the four 
elastic equations. Particulars of the bars of this truss are given in Table 2. 

TABLE 2 
Values of L and L/A for all members of the truss 


Lengtl 


Ba lil pane Panel I Panel II | Panel III | Panel IV 


12 30 } 20 15 30 
34 3 | 20 15 30 
13 50 5 60 100 30 
24 50 50 60 100 30 
23 40 48 40 24 20 
14 | 40 48 }O | 2 
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TABLE 3 


Elastic constants for panel groups 


Panels I and II Pane 


I it Be i Ld 

15 15 13 500 

{ I ) 15 13 500 

15 Oo I 250 000 

} 150 000 I 250 000 
3 64 24 35 400 
7 637 34°5 55 150 

| 463 JOO 620 550 


Structural analysis. Step 1. The results of the first step (determination 


of the values of S’ for the equivalent bars) are summarized in Fig. 8. 


For panel I alone, the values of L’(L)* are readily obtained (Table 3) from 
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ae: “ae 4 4 A 
7 y | Y | 
saiosiaiilh ‘ oer 
Fic. 10. ¢ aiaiia forces (step 3). ’ 
the given information. The sum, viz. 446,900, is the elastic constant C 


for the panel Hence, from (2 


( £46900 


D 279°3, 
L,)° (40)? 
nd from (7) 
D— Ls 18(279-3— 48) " 
S Leal =) 18(0-828) — 39-7. 
. D 279-3 
If this value is now used as the L,;, value in panel II (the other values 
eing as already given) a similar process will give for panels I and II 


combined the following values 


( 163.700 D 290 So $0(0-S862) 34°5. 

Similarly, the following are the values for panels I, II, and III combined: 
( 620,550 D SS Sos 24(0°938) 22:5. 

These values (of S,,) are thi ues quoted in Fig. 8a; those given in 
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Fig. 86 are obtained from the opposite end of the truss and this time they 


are the S;, values. 
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Fic. 12. 


Induced and total forces in bars (steps 5 and 6). (a) Forces induced from 
left to right. 


(6) Forees induced from right to left. 


TABLE 4 


Constants for equivalent panel III 


L’(L)* values jor 


conjugate jorces on 





B I L'(L}?? Left Right 
30 15 125 
30 15 13 500 13 500 13 500 
5 100 250 00 50 000 
5 100 250 000 250 000 
40 20°6 32 950 32 950 
40 +5 55 Iso 55 150 

Totals 615 100 318 650 296 450 


(= B) 
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ey Step 2. The constants for the four equivalent panels are shown in Fig. 9. 


They may be determined conveniently from a tabulation such as Table 4 


” 


















































+8 +14 14-2 
c < 
A 
r 
¥ ¥ \ 
' 
Fic. 12 (contd Induced and total forces in bars (steps 5 and 6). (c) Primary 
ym forces (collected from Fig. 11). (d) Total forces (a) +-(b)+(e). 


which deals with equivalent panel III and which is self-explanatory. It 
should be noted that both in Figs. 8 and 9 the dotted lines indicate 
members replaced by equivalent bars. 

Stress analysis. Step 3. The conjugate forces are obtained by statics 

» for the given loading and are shown in Fig. 10. 

Step 4. The primary forces in the members of each panel may now be 
written down on a sketch as in Fig. 11, using Tables 1(b) and 1(c), and 
finally collected together as in Fig. 12c. In Fig. 11, for panels II and ITI, 
two sets of forces are shown. These are due to the conjugate forces on the 
right and left respectively; they have to be added to give the total primary 
forces in the bars. 

Step 5. The forces induced from left to right are as shown in Fig. 12a. 


5092.9 . 
, E 
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The primary force + 9-62 (= P) in bar 23 of panel I induces a force in 
bar 14 of panel IT of amount P(—Z,,/D,) +- 9°62(—48/293-9) = —1-57, 
This gives rise to a set of internal forces in panel II, the value of the force 
in bar 23 being also —.1-57. In addition there is a primary force in this 
bar of +6-08 giving a net value of +-4-51 which induces — 0-463 in bar 23, 
panel III, and so on. The forces induced from right to left are given in 
Fig. 126. 

Step 6. The total forces in all bars are obtained by summation of 
primary forces (Fig. 12c) and induced forces (Figs. 12a and 126). They 
are shown in Fig. 12d. For chord members and diagonals this is a 
straightforward operation. For vertical members the primary forces are 
already included in Figs. 12a and 126, so that, for example, the total force 
in the middle vertical is obtained as 


1 4°51 —0-463-+- 1-064— 0-147 14-964, 


7. Conclusions 

1. The introduction of the concept of equivalent elastic systems pro- 
vides a ready means of solving pin-jointed redundant plane frameworks 
without the necessity of formulating and solving simultaneous equations. 
All computations can be performed by slide rule. 

2. The division of the analysis into two parts, of which one (the struc- 
tural analysis) is valid for all conditions of loading, is particularly valuable 
in cases where the effect of a number of different loading conditions has 
to be investigated. 
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THREE-DIMENSIONAL TURBULENCE AND 
EVAPORATION IN THE LOWER ATMOSPHERE, I 


By D. R. DAVIES (The University, Sheffield) 
Received 6 September 1949] 
SUMMARY 

The investigation is a first approach to the problem of turbulent diffusion from 

finite rectangular lake. A formal solution in terms of Mathieu functions is obtained 

r the problem of vapour distribution above and down-wind of the lake on the basis 

an assumed law of lateral diffusivity. The air-stream is assumed to be fully 

irbulent and directed parallel to the longitudinal length of the lake. The applica- 

tion of the formal solution to special cases awaits the complete tabulation of Mathieu 
functions. 


1, Introduction 

A THEORY to account for the rate of evaporation from plane saturated 
liquid surfaces into a fully turbulent air-stream was formulated by O. G. 
Sutton (1) and verified experimentally by Pasquill (2). The mathematical 
model of diffusion worked out by Sutton allows for the transfer of mass 
vertically only and, although quite satisfactory for assessing rates of 


evaporation, it does not give as good an agreement with experiment for 
the distribution of vapour above and down-wind from the evaporating 
surface. One improvement is suggested in the theoretical mechanism by 
the introduction of a Jateral diffusivity. This is discussed by the author 
na paper (3) in which an assumed simple power law of variation of lateral 
diffusivity with height is introduced to allow for this effect, and the 
resulting differential equations are solved to give the vapour distributions 
wer plane liquid surfaces which have parabolic boundaries: it was necessary 
to assume in this paper that the parabolic area extended indefinitely down- 
vind, so that it was not possible to obtain the vapour distribution down- 
ind of a finite area. 

The aim of this paper is to integrate these differential equations, subject 
to the conditions suitable for evaporation from a rectangular lake, of finite 
length and breadth, when the air-stream is directed parallel to the length 
of the lake. 

As in the problem of the parabolic area, a solution is first obtained based 
n the assumption that wind speed and diffusivity are constant. In this 
first approximation suitable transformations are introduced which reduce 
the fundamental diffusion equation to that of heat flow in two dimensions. 
This is solved by a method similar to that used by McLachlan (4) for 
the problem of heat conduction in elliptical cylinders. The method is then 


[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 
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generalized to the problem of diffusion into an air-stream in which the 
wind speed and diffusivity, both vertical and lateral, obey simple power 
laws. 

By an extension of another problem in heat conduction the vapour 
distribution down-wind of an evaporating area of any shape has been 
obtained. This completes the rectangle and also the previous parabola 
problem, on the basis of the hypothetical model of turbulence used by the 
present author. 


2. The equation of diffusion and boundary conditions 

If u denotes the mean wind velocity at height z, in the positive sense 
of the a-axis, and K,, K,, K, are the coefficients of diffusivity in x, y, and 
z directions respectively, the equation satisfied by the concentration, i.e, 
mass per unit volume, x(%,y,z) of vapour is 


xX 2X _ 2 | x x) 4-2 [x Xx], Ol gO], (2.1) 
Ct = ex xl Cyl * Oy| * dzl * a ; 


In this equation K,, K 


, and A, are such that the rates per unit area at 
which x is transferred across the planes x, y, z = constant, in the positive 
directions are —K,(0x/éx), —K,(éx/éy), and —K,(éx/éz). 

Steady conditions are assumed and K,, is neglected, so that the diffusion 
equation (2.1) becomes 

ies a] “ x22. (2.2) 
Cx Cy ~ CY Oz C2 

It is assumed that the vapour pressure of the evaporating liquid approaches 
its saturation value at the temperature of the surface as the surface is 
approached. This implies that the vapour concentration y assumes a 
constant value y, at the surface of the saturated area. 

The appropriate boundary conditions then are 


(i) x/x9 > 1, as e > O for values of x and y inside the area, 


(ii) lim{K.(éy/éz)} > 0 for x and y outside the area, 
0 


(iii) y> 0,2 > 0,z> 0, 
(iv) x>0, y>+0,2>0, 
(v) x>0,x> +0, z 0 


3. Wind speed and diffusivities constant with height above 
evaporating surface 


It is convenient to substitute non-dimensional variables defined by 


x = £aU, y = w(2K,), z= G/(«K.), (3.1) 





where 
consté 
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the where a is a scale factor of the dimensions of time, and U, K,, K, are 
wer constants. In these coordinates equation (2.2) becomes 
CY Cy a , 
RE an 4, (3.2) 
our o& on* 06° 
een 


with the boundary conditions 





ola ‘ ; 
the (i) x/Xo> LI > E> 0, oF Se, C0; 
(ii) €x/e¢ > 0, for € <1, yn? >c*, € > 0 and for € > I, all yn, > 0; 
(iii) y >0,€>+0,f> 0; 
(iv) x > 90, 0, ¢ 00, 
nse ‘ i 
nd where 2c is the breadth of the rectangular lake in (€, n, ) coordinates and 
i.e. | the length of the lake. 
Assume the existence of solutions of the type 
1 x = e~“*g(n, 6) (3.3) 
1) 
and substitute in (3.2) which becomes 
at e2 2° 
iv Ax sta |x: (3.4) 
ive Cc ~ Of os 
| This equation, together with the form of the boundary conditions, 
ion ae mee ‘ , 
suggests a transformation into elliptic coordinates (8,6) defined by 
nal ce cosh(B+-76). (3.5) 
.2) Assume xy = e“£B(B)H(8), (3.6) 
and substitute in equation (3.2) which is then separable into 
1es 
d*B 7 
mB nF Ac? cosh?8|B = 0 (3.7) 
a dp? 
' d*H — ‘ 
and ~-+|A—Ac* cos*0|H = 0, (3.8) 
dé? 
where A is the constant of separation. 
The boundary conditions become 
(i) x/x¥9 > 1, 1 € 0, 0 6<7,B->0; 
(11) 0x/C0 0. 1 é 0. O B 0. 6 QO, 7: 
ili) €y/e@ = 0 at all points on ground level for € > 1. 
Periodic solutions of Mathieu’s equation of ‘elliptic’ type (3.8) are 
ve given by 


ce,, (8, A) S Az” cos 278, 
r—0 ss 
defined and extensively tabulated by E. L. Ince (5); the symmetry about 


1) $= 47 cuts out the ce, of odd order. The solutions of the associated 
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‘hyperbolic’ type of Mathieu functions (3.7) have yet to be systematically 
tabulated. Definitions and references to existing tables are described by 


W. G. Bickley (6). The most suitable hyperbolic type of solution is given 
by W. G. Bickley in the above paper (equation 18.13) of the type 

2A Ee 1 oe — 

Sess 1)43" J,(ke-PV¥, (ke), 

A? - 

where k? = Ac?/4. Following the notation now in use these functions are 
denoted by Fey,,(8,A). They are known to converge uniformly in any 
finite region of the f-plane, including the origin; and they tend exponen- 
tially to zero for large 8 values. The existence of zeros of these functions 
when £ = 0 follows from their asymptotic expansions at large A values, 
Adopting the methods previously used by McLachlan (7) and Goldstein (8) 
in other physical problems, a formal solution expressed in terms of Mathieu 
functions is 


X — 1—Y¥ C,,ce,,(0,A)Fey,, (B, Ae, (3.9) 
Xo n=0 


is a function of A to be determined. 


2n 


where C, 


n 


At the surface of the lake 8 = 0, and when é > 0, y/x,—> 1. Hence 
Fey,,,(0,A) = 0. (3.10) 


Thus A has those values A,,,,, which makes Fey,,(0,A) vanish, i.e. the 
parametric zeros of the function. The function (3.9) clearly satisfies the 
condition that éy/e6 = 0 at @= 0, 7. At the up-wind edge of the lake, 
the physical conditions may be represented mathematically by making 
x/Xo > Oas € > +-(), 


Hence, > Dd Coy Feys, (B.A 


n=0 m=1 


)ce 


2n,m/ 


(9, Aon am) Ss (3.11) 


> 
2n 


This system of equations enables the C,,, to be evaluated. Use is made of 
the orthogonality theorem described by McLachlan (7, p. 176) and the 
development is similar to that used by him for the problem of heat 
conduction, except that integration with respect to B is taken from zero 
to infinity. A full explanation of this process is given in the general case 
described in section 4 of this paper. The formula obtained for C,, is 


| Fey,,(B, Aap )[ 242" cosh 28—A3"] dB 
a ; (3.12) 
| Fey3,,(B, Ao» m)[cosh 28—O,, ] dB 


0 


2n,m 


As A+ > At AP... (3.13) 


where 2) 


2n 
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The vapour distribution above a rectangular lake up to the second 


lateral edge is thus formally given by 


X 1 S > erent Fey,,(B, Aon m)C@an(9s Aon.am) X 
Xo gee * anpcey 
| Fey,,,(B, Asn. m)[ 242" cosh 28—A3"] dp 


i . (3.14) 
( Fey3,,(8, Aon m)[ cosh 28B—O,, | dB 


\ general theorem in heat conduction theory may now be used to give 
the vapour distribution down-wind of the lake. If v is the temperature at 


time ¢ at a point (a, y) the equation of heat conduction in two dimensions is 


ts Ke, nib (3.15) 
da? © by? 
where A is the heat diffusivity of the substance. 
By methods described by Carslaw (9) it may be shown that 
v,y,t) 
l 


He f(a a y' Ve j—y' AK p-((ev-2'P+Hy+y'} on} da'dy’ (3.16) 
7 


satisfies the conditions 

(i) v = f(x,y) at time ¢ Q, 

(ii) the solid is bounded by y = 0, this boundary being impervious to heat. 

The corresponding solutions for vapour distribution are immediately 
obtained by putting A l and t = (€—l). 

In order to calculate the distribution of vapour down-wind of the lake, 
the first step consists of evaluating numerically the distribution x = f(7, ¢) 
in the plane £ l, by using equation (3.14) above. Then for € > J, 


1 | f(’,¢ Te C-0')*}/4E-D 1 p—f(n—n’'P? 104-07 3/4E | dy'df’. 
(3.17) 
4. Wind speed and diffusivities variable with height above the 
surface: solutions up to the down-wind edge 
Following a previous paper by the author (3), a practical model of the 


diffusion mechanism in a turbulent air-stream is obtained by writing 


7 {+~\m 

(=) (4.1) 
U; =] 
K = «4,Urs, (4.2) 


K, - a, U} NM | (4.3) 
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In these expressions U, is the value of the mean wind speed at some 
standard height z,; a,, a,, m, and m are constants for a given turbulent 
state and may be worked out numerically from simple meteorological 
measurements. They are described in detail in the above-mentioned paper 
and in the papers by O. G. Sutton. The relationship between m and n is 
m = n/(2—n). A new development of the two-dimensional theory by 
K. L. Calder (10) yields different expressions for these constants and gives 
much improved agreement with observation, but the laws of variation of 
U and K, with height in his theory are similar to equations (4.1) and (4.2) 
above. The present work may therefore be easily adjusted formally to 
incorporate the new constants. 

The so-called conjugate power laws (4.1) and (4.2), although based on 
the assumption that the eddy shearing stresses are constant with height 
above ground level, are now generally accepted by meteorologists. As yet 
no other three-dimensional theory of atmospheric turbulence applicable to 
evaporation has been developed, and the form of equation (4.3) is chosen 
as it appears to be the only simple power law for K, which, when inserted 
in the transport equation, leads to a tractable partial differential equation. 
However, the quite significant improvement in the calculated rates of 
evaporation and in vapour distribution at the down-wind edge of a saturated 
parabolic area, quoted by the present author (3), gives physical support for 
the use of equation (4.3). Further experimental evidence to show that the 
‘z™ law’ may be used quite satisfactorily for evaporation problems is given 
in paper IT of this series. 

Equations (4.1), (4.2), and (4.3) are substituted in the transport equation 


(2.2) and the result simplified by writing 


& = 2, n = (UY /a, Z')*y, C = 2(Uf/a, Ze )yiemt4/(2m+1). (4.4) 
The fundamental equation becomes 
». OY 7) ne Ol. Ox ” 
hg ‘Xx : ue x 1 a xX ‘ (4.5) 
o& CC CC OnL ©7n. 








where s (2—n)/(2-+-7). 
Using elliptic coordinates »+7i¢ = c cosh(8+-7@), equation (4.5) becomes 
l 


9 9 





oe 
2 22 2 27 599 |X7 
c*(cosh?8 — cos?6) op? c@ } 
s . ‘ c CY 
+. P —— — : cosh B sin 6 — + sinh B cos 6 XY x 
c*(cosh?8 —cos*@) sinh f sin 6 cp 06 |’ § 
(4.6) 


This equation is separable by writing 


x = F(E)G(B)H(9). 
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Then it follows that F=e~%, (4.7) 
d? d e 
H L_sceoté H | (A —Ac? cos?6)H 0, (4.8) 
dé@2 dé 
d?G dG 
and => +s coth B — —(A—Ac? cosh*B)G = 0. (4.9) 


Equation (4.8) has been discussed by E. L. Ince (11). It may be shown 
that it is possible to expand solutions of this equation as a series of 
ordinary Mathieu functions. Ince begins with an equation in the form 

du 9(9—1) 


a+ 26 cos 2z 


_ a 0, (4.10) 

dz* sin-z 
which reduces to Mathieu functions when # = 0 or 1, and terms the 
periodic solutions as Associated Mathieu functions. He shows that (for 


y }) the solutions of (4.10) are also solutions of the integral equation, 


u(z) = p | e€ ©°S* 0°88 sin®z sin su(s) ds, (4.11) 
i 
the use of the symbol s in this equation being clearly distinguishable from 
its use in equations (4.8) and (4.9). 
3y using the orthogonality properties of the ordinary Mathieu functions 
it may be shown that 
pone > ay, C€s,,(2)Ces,, (8). (4,12) 
The functions ce,,(z) are known to satisfy the integral equation 


C0,,(Z) = pe ( e* CoS 088 Sgn. (sg) de. (4.13) 


0 


Substitute for e* ©°°* © 5, then 


if the ce, functions are normalized so that 
ce3.(z) dz = 2. 


] 
Hence Yo ‘ (4.14) 


: TfL, 
Ince denotes the even periodic Associated Mathieu functions by ce? (z), 
where the index n» has similar significance to that used in the ordinary 
Mathieu functions. Substitution of the expansion (4.12) in equation 


$11) o1VeS -_ 


ce,,.(z)ce,,(s) | sin’z sin”s ceg,,(s) ds. 
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The factors p, and p?, may be evaluated by methods similar to those 


given by McLachlan (7, p. 181). The equation for ce?,(z) may be put into 


the form 


e 


P : a . 3 9 . 
ce? (z) = Bwginez S ce,,(z) | ce? (s)ce,,(s)sin®s ds 


2n / 
7 << Hor . 
: 0 
ints FS baa. fs r 
sin’z > b3” ce,,(z), (4.15) 
= 
pe I i 
9 < off 
where i =n cey,,(s)sin’s ce,,(s) ds. 


T Mop . 
0 


This shows that it is possible to expand Associated Mathieu functions 
as a series of ordinary Mathieu functions, although it is not possible to 
obtain explicit expressions for the coefficients. In a numerical computa- 
tion use must be made of numerical methods on the lines followed by 
Ince (5). For the purposes of the present development the expansion 
given above is convenient. 

If w = u,(z)sin®z is substituted in equation (4.10), w,(z) is seen to satisfy 


1? l 
“41 1 99 cot 1 + (a—92-+ 20 cos 2z)u, = 0. (4.16) 
dz? - . az 

This is identical with equation (4.8)if 2% = s,a—9*— 20 = A,and 46 = —Ac®. 


Hence solutions of equation (4.8) may be written as a series of ordinary 


Mathieu functions H,,(0) = > b3" ce,, (0,2). (4.17) 


Since equation (4.9) may be obtained from (4.8) by substituting @ = if, 


then similar expressions must exist for equation (4.9). Write solutions of 


(4.9) in the form nwa 
Gan (B) = > cdf Fey.,(B, ). (4.18) 
t 


A formal solution to the problem is then 


X = 1— C4, (A)Hp,(0,A) Gy, (B, Ae, (4.19) 
Xo n 


where C;,,(A) is to be determined by the boundary conditions. 

It is now necessary to consider an orthogonality theorem for the 
Associated Mathieu functions on similar lines to that described by 
McLachlan. Let ¢,,(@,A) be a solution of the Associated Mathieu equation 


dg? 


]2 9.9 — 
“ape + | tn d,. cteos 29 — 27 |e = 19, (4.20) 


sin?@ 





and | 
soluti 


Int 


fa 


in 








hose 





: and let A,,,, be such that #,,(8,Anm) = 0 for B = 0, where #,(B,A) is a 
into ; 
solution of the corresponding equation 
dds o(P—1) 
a A,,,,, c2 cosh 2B — \ ub 0. (4.21) 
j mh2e |?” 
dp sinh?s | 
Write ¢ dh , then for (4,, »»;Anm) 
L.15 eC i , l l VW, 
— : “A (cosh =p COS 20) Ho L) . 9 ° > Ir m VU. 
ep c6 | sin?@_sinh?8 |J°”” 
(4.22) 
Similarly it follows tor (a A ) 
= a A I Ll hh, 
: =! ot | e2) (cosh 28 — cos 20) — (9 i ener > lr 0. 
ione ep ch | sin?@— sinh®B | | ~”” 
: to (4.23) 
ita- Now multiply (4.22) by ¢ 1.23) by f,,,,, and subtract, 
by C ( rae c ¢ 7) pe fc oo 
; >| > : | >p.r — >n,m 
s10n Cb Cp cp | co ob oo 7 
_ a c?(cosh 26 COs 20)(Anm a oe oe 0. (4.24) 
sty : . 
. Integrate (4.24) with respect to 6 from 0 to ~ and @ from 0 to 27, 
16 r ‘ ol ai 
)) cl : ‘ . hae 
| C ao a6 + | Coe et Sam SE | ABH 
: ‘ op cp |, 7 ~  o6 OF ie 
0 0 
\c?, - 
iry 2c7(A A.) | | (cosh 28B—cos 20)C,, m Spr IBdO = 0. (4.25) 
17) = ee 
lake d, — ce?(0) = sin®0 > b3” ce,,(0) 
. r 
4) ’ ) : % . 9 . f 
tp, and Y Fey? (8) = sinh’B > c3 Fey,,(8). 
ol t 
By using the asymptotic formulae for Fey,,(8) given by McLachlan 
5) (7, p. 220) it may be shown that Fey Y(B) >Oas B > OO. By virtue of this 
fact, the periodicity in 6 and the use of equation (4.28), the first two 
integrals vanish. Hence if p + n, 
9) | | (cosh 28—cos 26), mp» dBd? = 9. (4.26) 
0 0 
This holds also if p n, 1 m, but if p = n andr = m, 
1e€ ye : + 
“ { | [Fey®,,..(B) |*[ce8, (4) (cosh 28—cos 20) dBdé + 0. (4.27) 
n 
At the lake 


surface B 0 and when £ 


)) 














0, (x/xXo) :. 
S Feys,,(0, A) 
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Now apply the boundary conditions to equation (4.19). 
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Hence 
QO, 
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and thus as in section 3, A has those values A 
vanish. 

It may be shown that near the surface K, behaves like sin®*-1@ (3s > 1 
in the atmosphere) and that @y/éz is proportional to éy/é6. Therefore, for 
points at ground level outside the wet surface the function (4.19) satisfies 
the requirement that 


2n,m Which make each Fey,,, (0, A) 


: » OM; 
lim | K, £* 0 at 6=0,7. 
c>0, C2 

At the up-wind lateral edge of the lake the physical conditions give 
x/Xo > 9 as €— +0. Hence the system of equations to evaluate C,,, follow, 


1 = SY C,, Hy, (0)G,, (B). (4.29) 


n ™ 





Employing the orthogonality theorem by multiplying each side of 
(4.29) by 
H,,,(9, Aon, r) Gon (8, A, a ,)(cosh 28 — cos 20 )sinh29B sin299, 
and integrating with respect to 6 from 0 to 27 and with respect to 8 from 
0 to infinity, then the right-hand side vanishes except when p = n, r = m. 
Hence, 
l f Hi, (9, Aon m) Gan (Bs Aon m)[ Cosh 28 — cos 26|sinh?°B sin296 dBdé 


0 0 


0 Q7r 


[ [ | Fey, (B, Aon m) )/? ‘[ ce 3(0,A 


2n.m) "(cosh 28—cos 26) dBd@. (4.30) 
00 
The vapour distribution above a rectangular lake up to the second 
lateral edge is thus given in terms of Mathieu functions by 


Xx ee ee * at Taw 
l Zz > € — Zz bsy ce,,(6, Aen m) > C3 I ey o(B, Non m) ‘ ° 
Xo n=0 m=1 r t 


¥ > 
> br ce,,(8, a > Cc rt I e Mi 21(B; Aon m) ‘ 
: r t 


< (cosh 2 28 cos 26)sinh®B sin’@ dB d0 
((l> D3? €¢,(8, day,m) | Pe: c2? Fey,(B, A», am) ( (cosh 28—cos 20) dBdé 


(4.31) 
5. Wind speed and diffusivities variable with height above the 
surface: solution down-wind of the lake 
Using the notation of section 4, write 


ce= as (| Y, 


and x = Eroub(y, z), (5.1) 
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where the significance of the symbols y and z is clearly distinguishable 


from their use in earlier sections. Substituting these in (4.5) it becomes 


Cds Cus C 2s L é 
a l r. 1 al 
No W > Y 2 oe + « 


“\" ey Oz Oy? 2° oz 


=) : (5.2) 


Cz 


Putting ud Y(y)Z(z) in (5.2) and separating the variables, we have 
‘ SP 


r wr’ ” ((1/2*)(d/dz)(z8Z’)-4 p2Z'| 
r 2 . Z | 
where accents denote derivatives. It follows that 
Y"+dyY’ BY (5.3) 
und ! : (z8Z’)+-42Z' rZ, (5.4) 
where a—fB Ng. 
A suitable special solution of equation (5.3), when £ 4, is 
Y = e-v'it, (5.5) 
To fit the conditions of the problem, choose 7, 1, then 
x B No 1. 


Rewriting equacion (5.4) in the form 


Z" Z' 1(zZ'+Z), (5.6) 
ind substituting a series solution in the form 


a= >a.2, 


the relationship between the coefficients of z in the general case is found 


to be 
(n 1) — 
a "eee (5.7) 
») a” - 
2n(n ]+-s) 
The modulus of the ultimate ratio of convergence is z?/2n, so the series is 


absolutely convergent for all values of z. In the present problem even- 
power series solutions are required, and in the following work these will 
be denoted by Z. 

To obtain the vapour distribution down-wind of the lake the method 
described by Carslaw (9) may be followed. Let &, €—l; then suppose 
the distribution in the vertical plane at €, = 0 has been calculated from the 
solution in section 4 or is known from experiment. Suppose that at €, = 0, 
yx = f(m,¢). Using the special solutions of (5.3) and (5.4) found above, 


- 0 must be of the form 


wit > vat. ol [StS aye . 
ié Z\- — | eAn-7 1517 = _ dy dl - (5.8) 
Vey vey 


the vapour concentration at any point (7, ¢) for €, 


Arr 
x , | | st C )i¢ 


. 
00 


where A is a constant to be determined. 
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Consider the first term in the integrand. It is sufficient to show that this 
may be reduced to }f(n,¢) for € > 0 when €, > 0. Write 


n—n = —vey 
and C—C’ = —vé,2. 


The first term in (5.8) becomes 


A [f(a-+véry, C408, 2)e¥*8Z(2) dyde. 
“oc 0 


Let €, > 0, then this reduces to 
| Zdz = Af(n, 0)2n*F(s), 
0 


Af(n, f) ( e-v'lt dy 


x 


where F(s) | Z(z,8) dz. Hence the final solution for € > 1 is 





l Pact oe 2 
mn’. —(n-'F/4E—-D) ZF} _> = al 
4a F(s)(€ =a | Sy c| J 


LW(Eé—Z)J 
Lon-nving-ng {| S15" H In'dt’. (5.9 
wwe: 


6. Discussion 

It must be borne in mind that the present paper describes a mathematical 
solution to the problem of turbulent transfer from a rectangular lake, on 
the basis of an assumed power law of variation of lateral diffusivity with 
height above ground. The application of this and other assumed laws to 
diffusion from a continuous point source and to evaporation problems is 
discussed in paper IT of the series. 

[t has been shown experimentally by Pasquill (2) that Sutton’s theory 
is adequate to account for the rate of transmission of heat from a flat 
surface by turbulent action. It therefore appears probable, if the tem- 
perature difference between the surface and the air-stream is not large, 
that the mathematical solution obtained in this paper may be applied to 
give the distribution of temperature above and down-wind of a finite 
heated rectangular plate kept at constant uniform temperature, sur- 
rounded by a plane surface impervious to heat. In this connexion y in the 
evaporation problem should be replaced by pC,,t, where C, is the specific 
heat of air at constant pressure, ¢ is the difference in temperature between 
plate surface and free air-stream, p the air density. 


The work of computing numerically the general solution giving the 
vapour distribution above a rectangular lake cannot be undertaken until 
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the Fey,,,(8, A) functions have been systematically tabulated. The complex 
nature of equation (4.31) shows that even when these are available much 
labour will be involved in the evaluation, so for the time being the solu- 
tion, like many other problems involving Mathieu functions, has a formal 
value only. 

However, the analysis of section 5 is applicable to any evaporating area 
and gives a method of completing the problem of vapour distribution from 
, parabolic area (3). Here the solution involves functions which may be 
readily computed; and the vapour distribution, in a vertical plane perpen- 
licular to the axis of the parabola at the down-wind edge worked out 
numerically. The formulae of section 5 now apply and the full effects of 
the assumed law of variation of lateral diffusion on the calculated distribu- 
tion down-wind of the area may be estimated. It is proposed to carry out 
this numerical investigation as time permits. 
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SUMMARY 


This paper gives a determination of a power law for the variation of the coefficient 
of lateral diffusivity with height above ground, and contains a discussion of the 
application of a previously used law to evaporation problems. The index in the 
power law is determined by comparison of theory with well-known results giving 
the variation of peak concentration (i.e. weight of suspended matter per unit volume) 
with distance down-wind from a continuous point source, such as a single static 
smoke generator, producing a cloud of air-borne particles of matter. The applications 
of this empirical law of diffusivity, and of one previously used, to evaporation 
problems are compared and discussed. 


1. Introduction 


THE problems of evaporation and diffusion of vapour, into the earth’s 
turbulent air-stream, above and down-wind from saturated liquid surfaces 
of parabolic and rectangular shape have been solved by the author. A solu- 
tion of the former problem was published by him in 1947 (1) and a solution 
of the latter appears as paper I of this series (2). The two-dimensional model 
of turbulence used by O. G. Sutton (3), and recently perfected by K. L. 
Calder (4) was adopted and a three-dimensional mechanism set up by the 
introduction of an assumed law of lateral eddy diffusivity (see 1, equation 
4.4). No other three-dimensional theory, applicable to evaporation prob- 
lems, exists, and the choice of power law for the lateral effect was governed 
by mathematical convenience and physical plausibility; it was supported 
by the very significant improvement (compared with two-dimensional 
theory) in the agreement between theory and experiment for rates of 
evaporation and diffusion of vapour (1). 

For the purposes, however, of obtaining a solution to the problem of 
diffusion of air-borne particles of matter from a continuous point source, 
such as a single static smoke generator, knowledge of a more exact power 
law is necessary. At the present time it appears that the only method of 
searching for this law must involve the following procedure: (i) one must 
accept the experimentally established laws of variation for wind speed and 
vertical diffusivity from two-dimensional theory; (ii) one must assume 
possible laws of variation of lateral diffusivity with height; (iii) one must 


{Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 
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obtain the corresponding continuous point source solutions and (iv) com- 
pare with known results of experiment. 
This method has been followed in this paper and an appropriate law of 


variation has been determined 


2. The equation of diffusion, incorporating the lateral effect 

If denotes the mean wind speed at a vertical height z above ground 
level in the positive sense of the x-axis of rectangular axes Ox, Oy, Oz, and 
K,, K. are the usual coefficients of eddy diffusivity in the y and z directions, 
the equation satisfied by the concentration x(7,y,z), ie. the weight of 


re liquid vapour per unit volume, is 
dsp yx? lK x % cea (2.1) 
5 Ca cy| “oy Cz "2 | 

ime 
tati vhen steady conditions are assumed and diffusion in the 2 direction is 
lons neglected. 
\t this stage use is made of the so-called ‘conjugate power law’ which 

yppears to be generally accepted by meteorologists, namely 

( 2 \m 
7 =(- (2.2) 

th’s i “a 
1Ces ind Kk a, Ur, (2.3) 
dlu- In these expressions U, is the value of the mean wind speed at some 
ion standard height z,; a., m, and n are constants for a given turbulent state 
del and may be evaluated numerically from meteorological measurements. 
LL. They are described in detail by O. G. Sutton (3). The relationship between 
the mand nism = n/(2—n). The recent development of the two-dimensional 
ion theory by K. L. Calder (4) yields new expressions for these constants and 
ob- gives much improved agreement with experiment, but the laws of variation 
ned of U and K, with height in his theory are similar to equations (2.2) and 
ted 2.3) above. The present work, therefore, incorporates the new constants. 
nal As equations (2.2) and (2.3) represent a mathematical model of diffusion 
of which agrees well with experiment for two-dimensional problems, it appears 

that a three-dimensional model will best be defined by accepting (2.2) and 
of (2.3) and writing RK a, Ul-nze, (2.4) 
ce, 
ver where « is an index which has to be determined. 
of The coefficient a, will depend on the index a, as K, has to be kept 
ist dimensionally correct. The expression for a, may be defined by analogy 
nd with a, as previously suggested by the author (1, equation (4.5)),and support 
me for this procedure comes from the fact that in the case of evaporation from 


parabolic areas it yields results in accord with experiment. 
5092.9 F 
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Substitution of equations (2.2)-(2.4) into (2.1) leads to an equation of 
the form ( 


gmoX _ 4 & at-nX) ry (e952), (2.5 


Ox 02 02 oy\ cy 


we 


where A = a,z}"/U? and B = a,2?"/U}. 


3. Evaluation of an empirical power law of variation for K,, with 
height above ground, from well-known results for a continuous 
point source 


ne y 28 
Write —- = t, —__ =v (3.1) 
gla ala 
and xX = xraf(u, v), (3.2) 
where the constants f, g, and my are to be determined, and the equation 
(2.5) becomes 


ou Ov ou" 


ano-l ous Cus , o.. On 
Ny coh —- 9 (m ¥ +) = B(vala)(a—m)By(no-2la) am 


‘i 
Tee (ics @)(28—2m- 1)/By(no—2 a) : “L. 


dv? 

; . Ous 
+ B(p— m)(valla)(B-2m—1)/By-(no—1/a) ae ‘ (3.3) 

Cc v1 
The partial differential equation (2.5) involving three variables x, y, and 
z may now be reduced to one containing two only, w and v, by arranging 
the indices of x on both sides of equation (3.3) to be equal to (n)»—1) and 

thus cancelling out 2, i.e. 


x—m 2 
Bq +N)— - = %—l, (3.4) 
2 2m—1 2 or 
p Bq No—- N —1, (3.5) 
— 9m — 
and p—2m ee —-=n—1 (3.6) 
Bq 0 0 
From (3.4) and (3.5) we find that 
B = (at+tm-+1)/2, (3.7) 
and gq = 2(2m+1)/(a+m-+1), (3.8) 


and these values also satisfy equation (3.6). 


Suppose now thatthere is a continuous point source, at the origin of coordi- 
nates, emitting @ grammes per second of air-borne particles of matter. 
If the strength of the source is constant, then the flux of matter across 
an infinite vertical plane must be equal to Q at all distances down-wind. 
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Hence yx must satisfy the following conditions: 


(i) ( ( u(z)y(x, y,z) dydz = Q for all values of 2; (3.9) 
20 

(ii) x >, E>n>C>0; x > 0, E€>0, »>90, ¢(>0. 

(3.10) 


Now x is a function of the form 
x = Caraf(u, v), 
where C is a constant, determined by equation (3.9). By using (2.2) and 
(3.1), equation (3.9) becomes in (wu, v) coordinates, 
CU, 


Q~m 


P24 


1 gino +1/a-+Hm+1)/Pa b(u, v)v+1-PB dudv = Q. (3.11) 


The coefficient of x is to be zero, and hence 


No [1+-(m+1)/B]/q. (3.12) 
Substitution of the values of 8 and qg given in equations (3.7) and (3.8) leads 
to No (a+3m-+-3)/(4m-+-2), 
or \ 2(2m-+-1)n9>—3m—3. (3.13) 


At this stage, if a satisfactory experimental law was known giving the 
variation of peak concentration (i.e. along the x-axis) with distance from 
, continuous point source for various degrees of atmospheric turbulence, 
such a law would enable n, to be determined in terms of m, and « follows 
from equation (3.13). Unfortunately, as indicated by O. G. Sutton (5), 
reliable data on diffusion over the whole range of atmospheric stabilities 
is not available. However, as mentioned by O. G. Sutton in this paper, 
it has been established experimentally, in adiabatic conditions, that the 
central or peak concentration from a continuous point source decreases 
vith distance down-wind (x) according to the law x < x~?'6, i.e. my = — 1-76. 

Substituting this particular value for n,) into equation (3.13) and taking 
the value of m as } leads to « 1-09, in the particular case of adiabatic 
conditions; hence in these conditions K, appears to increase with height 
ta rate slightly greater than the linear law. 

Although it is not possible to obtain an exact solution to the problem 
f the determination of « to cover the whole range of atmospheric stabili- 
ties, an approximate theory may be developed. O. G. Sutton (5) gives a 

ntinuous point-source formula based on his ‘statistical’ theory and from 
this work the law of variation of peak concentration is written 2-@-™, 
The index has been fully verified in the neighbourhood of adiabatic con- 
ditions. As indicated by Sutton, this theory is based on the assumption 
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that wind-speed variation with height may be neglected and this leads to 
only a small error in practice. However, if comparison is to be made with 
the law x-@-”, the basic equation (2.1) must be treated with m 0 in the 
term U(éx/éx), i.e. the assumption must be made that wind-speed variation 
with height may be neglected. 


We therefore begin by writing equation (2.5) in the new form 


CY C OX c CX 
X— 4“ (zgi-mX) 1 BS (208X), (3.14) 
Cx Cz Cz | cy cy 

The treatment proceeds precisely as detailed above, but from this new 


equation. The value of f is found to be given again by equation (3.7), ie. 
B (at+m-+1)/2, 
and the new expression for q is 
q 2(m+-1)/(a+m+1). (3.15) 
When the variation of wind speed with height is neglected, equation 
(3.12) for n, becomes Me (1-+1/B)/¢. 


Substitution of equations (3.7) and (3.15) into this formula for , leads to 


Ny = (a+m-+3)/2(m-+1), 
or Y 2(m+1)ny,—m—3. (3.16) 
The relation 7, (2—n) is given by O. G. Sutton’s statistical theory, 
and substitution in equation (3.16) gives a 1—m. On the basis of this 


approximate theory the empirical power law of variation of A, with height 


therefore appears to be K a. Ui-nei-m 
y 7] 1 7 ? 


which is the same law as that used for K,. 

In the actual flow of air over a rough surface the physical conditions are, 
of course, very complex and probably different laws of variation of wind 
speed and diffusivities hold in different height ranges, especially near the 
ground. The best that we can hope to achieve with a mathematical theory 
at the moment is to assume that our model of diffusion is true over the 
whole height range where diffusion is expected to take place down to z = 0. 
Complete agreement with experiment justifies this step in two-dimensional 
problems; in three dimensions agreement is necessarily obtained for peak 
concentration variation, but a complete test awaits further theoretical 
work on cloud widths from a continuous point source. When the appro- 
priate values of a, g, and m, are substituted into equation (3.3), it leads to 
an equation, in w, v coordinates, which is not separable, and it appears that 
the full solution for a point source requires a numerical evaluation of the 


% function. It is only in the particular case of « = m that the equation is 
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separable and the continuous point source solution may be obtained in an 
explicit mathematical form, as given in Appendix I. The application of 
the case «a m to evaporation problems is discussed in the following 
section. 

4. The evaporation problem: use of K, = a, Uj-"2" 

[t is clear that the two-dimensional evaporation theory worked out by 
0. G. Sutton (3) and K. L. Calder (4) is only applicable if the cross-wind 
width of the evaporating area is large, except at points on the central axis 
of the area, not too far down-wind. At points just outside the area con- 
tained between the lines running directly down-wind, from the longitudinal 
edges of a rectangular evaporating area, for example, the two-dimensional 
theory obviously gives zero concentration of vapour, whereas it is well 
known from experiment that a considerable concentration is in fact pro- 
duced. This is due to the fact that near the surface the lateral eddies have 
a greatel energy than the vertical components. 

The author has succeeded in solving the problem of evaporation and 
vapour distribution over and down-wind of areas of parabolic (1) 
and rectangular shape (2), using a model of turbulence given by 
K a, Ut-"z™, The parabolic solution was expressed in terms of easily 
computed functions and the rectangular solution in terms of Mathieu 
functions which are not yet completely tabulated. From the table in 
\ppendix 2 it is seen that for a continuous point source solution the 
experimental variation of peak concentration down wind follows the law 

176 whereas in the case A, proportional to 2” it is x-?°; both in adiabatic 
conditions. For large values of x the error introduced by using the 2” 
law is clearly very considerable. However, an evaporating area may for 
mathematical convenience be considered to be equivalent to a very close 
grid of point sources, and at a point over the area the major portion of the 
concentration is due to those sources near the point under consideration 
and hence the previous solution for the vapour distribution over parabolic 
and rectangular areas may be expected to lie very near the truth. This 
also applies to distances not too far down-wind: the actual limit may only 
be determined by experiment. In support of this expectation the great 
improvement in the theoretical vapour distribution over a parabolic area, 
indicated by the author (1), must be quoted. In addition the results of 


experiments in the field, arranged by the Ministry of Supply, over the 


down-wind edge of a contaminated parabolic strip of short grassland, 
10 yards long, indicate as closely as can be expected in field trials of this 
type that the ‘2’ law’ produces a fully sufficient degree of lateral diffusion: 
this effect was particularly noticeable at the outside boundary of the area, i.e. 
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at a distance of 5 yards from the axis of the area. The results of a typical 
experiment are given in Table 1. The method used for the theoretical] 
computation has been described previously by the author (1). Direct 
comparison between theory and experiment is not possible since theory 
gives the concentration y(, y,z) at any point above the area in terms of Xo» 
the ground concentration on the evaporating area. As y, is not known, 


TABLE I (a) 





Lateral | 


dict > 
stance } 
2 Height of sampling point 


from the Type of 
axis result et = | 24” 18” 
o yd. Experimental 1700 «6| O35 «|| «0°06 Oro! 
Theoretical, based on K,,&2"| 1-00 | 0°38 | or 0:02 
Theoretical, based on K,, = 0} 1:00 0°46 016 O04 
| | 
25 yd. Experimental 068 | 0:26 | o04 | 0004 
Theoretical, based on K,, 2" | o-95 0°33 008 =| 0007 
| Theoretical, based on K, ©} O95 | 0°39 | orl 0°02 
5°0 yd. Experimental O15 | 005 | ool 0'003 
Theoretical, based on K,, C2” | 0°37. | o-og | o-02 | 0-003 
Theoretical, based on K,, = 0| 0 ° lo | 0 
| 
7°5 yd. Experimental Oro! 0004 | 0003 | 0002 
Theoretical, based on K,, &2""| 0-003 | 0-002 | o-oor | o 
Theoretical, based on K,, oo] o ° ° 1 o 





Theoretical and experimental results giving mean values of relative concentrations at 
various points in a vertical plane above the down-wind edge of a parabolic evaporating area, 
initially contaminated with aniline to a density of 10 g.m.~? The area was 10 yards in length 
and 10 yards in width across the down-wind edge, with axis directed along mean air-stream, 
and vertex pointing up-stream. 


TABLE I (b) 


Observed meteorological data, measured over the period of the experiment. 
comparison is made by expressing both theoretical and observed con- 


centrations in terms of the concentration at a height of 1 in. at the centre 
of the down-wind edge of the area. 


If the vapour distribution needs to be calculated at large distances down 
wind of the area, then the solution based on the law making K,, proportional 
to 2” will enable the distribution to be evaluated across the down-wind 
edge of the area—the first necessary step—and the method indicated in 
section 5 of the author’s paper on the rectangle problem then adapted for 
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the more exact K, law, using the substitutions (3.1) and (3.2) of the present 
paper. The problem, however, of evaluating numerically the function ¢(u, v) 


will first have to be solved. 


5. Discussion 

The dependence of the theoretical variation of peak concentration from 
a continuous point source with distance has been worked out; and an 
empirical law of variation, with height, for the lateral coefficient of eddy 
diffusivity chosen to give agreement with well-known results. This new 
law differs from the law previously assumed (for mathematical con- 
venience) in three-dimensional evaporation and diffusion problems. It 
appears that the previous law yields a point source solution providing 
insufficient diffusion, in the direction of the mean wind, at large distances 
from the source; but at small distances the error is small. Experimental 
results are quoted which indicate that the previous law gives good agree- 
ment for vapour diffusion at heights over an evaporating surface. This 
agreement is likely to continue until the diffusion at large distances down- 
wind of the area is to be considered. For the solution to this particular 
problem and for a full solution to the continuous point source problem, 
the function ys, defined in the present paper as a particular solution of a 
second-order partial differential equation in two variables, needs to be 
evaluated numerically. 

[ wish to express my thanks to Professor L. Rosenhead, F.R.S., for his 
interest in and criticism of the work. Acknowledgement is made to the 
Chief Scientist, Ministry of Supply, for permission to publish the contents 
of Table 1. 
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APPENDIX I 
Continuous point source solution in the case K, proportional to 2™ 


The equation (2.5) may be simplified by means of the following transformations : 
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The fundamental equation becomes 


ot 2 [per], 2 [ea 
~ €€ cn cn cGk” cg 


wu 


where s (2—n)/(2+7n). 
n C 
Now write ‘ U, = v, 
VE vé 
and x ENG ( a, 2), 


Substitution of (3) into (2) leads to the equation 


/ j > , ‘ 
. l Cus cw C7 1 ¢ Cus 
Naw { u b) } { a - 
o¢ > > . . 
ys Cu Cv. cu v" CV cv 


To solve, put us Y(u)Z(v), then 


ea , | Z J 
Hence it follows that AuY’ BY, 
Ld 7 ai 
and ; 7, (vZ’) luZ vZ, 
where x B Tn. 


In the above ie; n: C) coordinates, the condition (3.10) becomes 


K | | xgvemdndt = Q 


‘ ‘ 2 U, (a,,2 2 ] l /a.2™ 14+1)/272/(27 
wher K iz ee [ ol ) (4:2) 


(2m 


A suitable solution of (6) is 


when 8 j 


The equation for the Z function then becomes 


l d 
3g, (S’) lvZ (4+n,)Z. 


Using equation (3), « quation (8) becomes 


ae ae oe ‘ 
KE" | | p(—h, $.) cuem+) dpat = Q 
ee ee NS NE ‘ 
4 U0 


Now use (3) in (11) and it leads to 
KE ot lt1/(2 1 | | uk(u.v)el 2 l) dudwv QO. 


Hence ) (}-- de). 


where s, (2m 1). Equation (10) becomes 


This equation has the particular solution 


Z e—v/4 eS /4€, 


(6) 


(9) 


(10) 


(11) 


(12) 


(13) 
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The solution in (€,7,¢) coordinates is thus of the form 
ee 
9 x ; xpi — (n° + ¢*)/45}5 (14) 
re A is a constant easily determined from (8). Equation (14) satisfies conditions 
10 
It may be seen that this solution for the case a = m may be arrived at immediately 
ssuming a trial solution of the form 
{ . 
\ j xp (ay? ra VE} Be (15) ' 
ting in equation (2.5 | comparing the indices of separate terms, the 
ropriate valu I > 7, & ma I determined. However, although the method 
1 abo rather involved comparison with the trial method, it may 
the ger lines along v the more difficult solution may proceed for 
when A taken to be ] portional to z! 


APPENDIX II 
Peak Concentration Variations (from Continuous Point Source) with Dis- 
. tance Down-wind, using Different Theoretical Lateral Diffusivities: Adiabatic 


Conditions assumed 


Va f | 
| Index 
K K U No 
(8 constant nstant constant | 1-00 


1-07 
1-40 
1-50 
zm 1-67 
ZN 1-75 


NOTE ON THE CIRCULATORY FLOW ABOUT A 
CIRCULAR CYLINDER THROUGH WHICH THE 
NORMAL VELOCITY IS LARGE 
By B. THWAITES (Imperial College, London) 
[Received 3 May 1949; revised 23 June 1949] 


SUMMARY 

In this paper we consider the circulatory viscous flow about a circular cylinder 
of infinite length, the porous surface of which allows a normal velocity of fluid which, 
at any instant, is constant at all points of the surface. Asymptotic series in inversi 
powers of a suction Reynolds number are obtained, in steady and unsteady flows, 
for the distributions of circulation about circles concentric with the cylinder, pressure, 
and torque on the cylinder. Various comparisons are made with conventional 
boundary-layer theory. 


1. Introduction 
THE properties of circulation about two-dimensional aerofoils and other 
bodies have received a certain amount of attention recently. In particular, 
it has been suggested (1) that it is possible to set up circulation about any 
suitably rounded surface in a uniform stream and this has been achieved 
experimentally in the case of a porous circular cylinder fitted with a flap. 
Whether or not, on removal of the flap, the circulation at infinity remains 
constant or decreases is important, and experiments are being designed to 
investigate the matter. These experiments will include attempts to obtain 
circulation about a fixed cylinder in otherwise still air. 

The flow considered in this paper is the radial symmetrical flow outside 


a fixed porous circular cylinder on which is superposed a distribution of 


circulation. Exact solutions have been given for such a steady flow (2) 
which show, inter alia, that the suction velocity must exceed a certain 
minimum value for a non-zero value of circulation at infinity: for otherwise 
the vorticity will diffuse sufficiently rapidly through the fluid to destroy 
this circulation. No general solution in finite terms exists for unsteady flow. 

The principal object of this paper is to find general solutions of flow in 
the form of asymptotic series in a parameter involving the suction velocity 
from which the velocity and pressure fields and torque on the cylinder can 
be calculated in cases where no general solution exists. These solutions not 
only yield interesting analogies to conventional boundary-layer theory (3) 
but also will provide direct comparisons with experiment. It is shown that, 
for values of the suction velocity appropriate to this asymptotic theory, 
the circulation at infinity remains constant. 


[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 





NOTE 


2. No 
Cyli 
velocit 


identi 


where 
If t 


value 


in wh 


boun 
the a 
of mi 


whet 
distr 
span 


Th 


hav 


























NOTE ON CIRCULATORY FLOW ABOUT A CIRCULAR CYLINDER 75 





2. Notation and equations 
Cylindrical polar coordinates (r,¢,z) are used, with the corresponding 
velocity components (v,u,w). If derivatives with respect to ¢ and z are 


identically zero and w = 0, the equation of continuity gives — (vr) = 0 or 
or 
v¥ = constant 1-94, Vo = Uo(t) being the normal velocity at the cylinder 


a. The d-equation of motion is 


ou cu. we ( 24 1 du “) 
a) v +-———-—], 


ot or r Or? ror r= 
nd . P , ‘ : 
hick where v is the kinematic viscosity. 
ersi If the circulation k = 27ur is taken as the dependent variable and the 
ows, value given above for v is used, this equation becomes 
sure, 
ek ck rok 
ona r—_.t(R—1)——-—=0 (1) 
or Or vdt 
. i : Vy , . : 
in which R is the Reynolds number - If the cylinder is at rest, the 
ie 
her atl r 
lar boundary conditions are k = 0onr = a, andk > k,(t) as r > 0, k,(t) being 
a the arbitrary time-distribution of circulation at infinity. The r-equation 
of motion reduces to 

ved 
; oR ie r3 Op 
lap. vr- v5 ar = + — A 9 (2) 
‘ins ot 4a p or 
l to where p is the density and p is the pressure of the fluid. This gives the 
ain distribution of pressure. Also p— py as roo. The torque M per unit 

span exerted by the fluid on the cylinder is given simply by 
ide 

ok 
of M pal) (3) 
r=a 
(2) ~— ; 
' where pu is the coefficient of viscosity. 
ain : 
‘a0 3. Steady motion 
‘OV (i) When the flow is independent of time, ¢, solutions of (1), (2), and 
WW. 3) are given in (2) and are 
in \ R-2 
kk = kp! 

Ib) | r| { 


l ke p\/a 2 1 /a R-2 l a 2(R-—2) 4n2v, a2 ; (4) 
ae emer ee 
om” \ ae r R\r R—1\r (2 





(2) 
at. M hey pel R 2) 
ys These solutions show that for conditions at infinity to be satisfied we must 


have R > 2, or that steady motion is impossible when —v, < (2v/a). They 
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also show that the circulation tends to its value k, at infinity slowly, and 
in particular more slowly than the exponential manner in which the 
velocity in a boundary-layer on a plane surface tends to its stream velocity, 
However, it can be seen that the circulation increases rapidly near the 
surface of the cylinder, the rate of increase increasing with R. But these 
solutions, as they stand, give no indication of particular behaviour when 
R is very large or of the limit of the circulation distribution as R—->o, 
However, these conditions may be approximated to in experiment and it 
is the purpose of this paper to display the behaviour of solutions in 
asymptotic series in R. 


(ii) Now the circulation-loss layer, or the layer near the cylinder’s 
surface in which the circulation increases rapidly to near its asymptotic 
value, varies with R and some assumption must be made about the 
variation. The vorticity which is diffused outwards from the surface is 
convected inward by the suction velocity, and so when the suction 
velocity is large, the total distance through which vorticity is transported 
outward from the cylinder is inversely proportional to the suction velocity, 
Thus, on the hypothesis that the circulation-loss layer is inversely propor- 
tional to the suction velocity, we put 


rT a- 


na/R, (5) 


» being a non-dimensional measure of radial distance y from the surface 
of the cylinder. Then » = (—vyy)/v and is a familiar parameter of con- 
ventional boundary-layer suction theory. If this form (5) for 7 is inserted 
into the solutions (4), it is possible to expand each expression in a series 
in inverse powers of R by elementary means to obtain: 


k keh l e-7) — “ (y 1 4)e-7 _ oT (37? | 16% 1 24)e-7 — 
4 ] 
7) 9 _ | 
2164 4)e-7-LO 
48 R8 (7 dic silil (7) 
(Po—p)r* 1 ko */ I " (4e-7—e-*”) (6) 
Lpvra? 2rv} | R? R8 


(24°+8n)e—(n?+-47-44 nie 4 
R4 R°}) 


- pak, V9 R 


These are the required asymptotic expressions. The first series shows that 
as R00, 


the distribution of circulation (in circles concentric with the 
cylinder) is precisely similar to the distribution u/V = 1 


evoul” which 
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and gives the velocity wu in a boundary-layer at a point where the stream 
the velocity is U and the suction velocity vy is very large (3) and negative. 
“ity. The second series shows that the pressure differences in the field are 
the proportional to v2 and the last that the torque is proportional to vo, a 
hese result which follows, of course, from the assumption that the circulation- 
hen loss layer thickness is inversely proportional to vp. 

Pin iii) The series in (6) can also be derived from the appropriate equations 
ys vith » as an independent variable in place of r. For example, in the 


steady motion case under consideration, (1) becomes 


er’s ( 7) te (12) 0. 
\ R = R CY 


tic 
the und on the assumption that 
is I bs 
te k = hol foln) +5 fil) + pp fol) to) | oi 
ted 
f,(0) f (0) J tap) (), n l, fi (2) J 
by 
ie successive equations are obtained for f,(7), f;(7) and so on, which can be 
ilved to give the series already obtained in (6). This procedure therefore 
; merely checks the simpler method of obtaining the asymptotic series from 
"7 the exact results. However, in any case where an exact explicit solution 
ce is not known, this procedure is the only possible one, and it is now applied 
n- to the case of unsteady flow. 
“(l 
- 4. Unsteady motion 
(i) If the circulation, ky, at infinity is not constant, then it can be easily 
leduced from (1) by considerations of the behaviour of k as r > oo that R 
must be less than 2; and when R < 2 it appears that the distribution of k, 
vith ¢ can be arbitrarily assigned.+ This presents a difficulty in the 
determination of the decay of /, when initial conditions are postulated. 
However, an asymptotic investigation for large R is unlikely to be valid 
for R < 2 when in any case the transformation (5) to 7 is inappropriate. 
) The case R 2 only is therefore considered when /, must remain constant 


and R or the suction velocity varies with time. It may be remarked here 
that, from (1), the conditions at time ¢, depend only upon what has 


ck Ck ine 4 ° ° 
With R l 1 particular exampk 1) becomes :v—,. This equation occurs in 
ct or 
etheory of heat conduction and also is of the form of the so-called outer solution of Karman 
and Millikan’s method of dealing with boundary-layer flow on a plane boundary. The equation 


in be solved in terms of the variable and in particular, a series solution is obtainable 
Dt 
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happened for ¢ < ¢,: this corresponds to the conventional boundary-layer 
theory result that conditions at any point depend only on what has 
happened upstream of that point. 


(ii) When R > 2, equation (1) possesses no general solutions of the 
simplicity of the steady solutions of § 3; investigation of the behaviour of 
the flow when R is large (and not constant) is not as simple as in § 3, 
therefore, and must proceed by the use of series such as (7). The asymptotic 
transformation 
R= R(T) = ae i = = T, o1+3), S(T) = dl (8) 


applied to equation (1) gives 


o*k 1\ ok l ok 
“wee the PEE eco 2) S—|=0. 9 
( , 1) + ae ( TR at (9) 
If the series (7) is inserted into (9), the f’s now lies ne of » and ¢, 
the following equations are obtained: 


Soy ofo oe 





— 
7” on 
Of, a, fy, 
On? on 1 Gy? on 
fe , fe _—— fh ch, ¢ Y rt eR 
on*® On 19,8 ont oT : 
re oe > (10) 
Bit t= 1 at Si ot! 92 —8f,ty bey th 
nn oN on TT 
OSs Ss = aie et + Sy och r 
On“ On én” on én 
C Ofs 
+ 8 22_ Snf,—28f,-+9: =o 
on oT 


and so on. These equations can be solved simply by separation of vari- 

ables; fy and f, are functions of 7 only, and are given in the first two 
terms of the series for k in (6); f, and f, can be written as 

1 Y p 4 

f(y, T) = J2(n)+ Shg(n), fs(n, T) = 93(n)+Sh4(n), (11) 

the g’s and h’s being functions of » alone. For example, the equation for 

h, is - , = 

i eh, , Oh, Cf 

——= — = 7 9 


ins BO ha(0) = hg(co) = 0. (12) 
C n~ C 1) Cc 1) 


The particular integral of f, clearly involves S, S’, and S?, and so we put 


fa(n, T) = Ja(q)+ Sha(n)+ 87) 4(n)+S’kq(n), (13) 
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aver andf,, m > 4, can be expressed similarly in terms of combinations of S 
ha and its derivatives. The g’s are given already in (6) (the steady case in 

which S Q). 

To the series in (6) for k must be added the following terms: 
} 5 

the 
ee kn(n+2)e78 1 9(3y°+327?+96n+192)e78 O 1 (14) 
§ 9 hk 24 R R)’ 
So 


otic ind so, from (3), to the series for //(— pak, vy) the terms 


S 8S l 
== +0 } (15) 


re Ve 
(8) |, ' — — 2 ; -™ 
[he effect of variation of R on pressure p is found from (2). First there 
is the addition to p due to the term vr? ‘ . which is easily calculable and 
c 
(9) means the replacement of Po . in (6) by 
Plea 
. (Po— P+ (pv*/a*) SR log(r/a)}r° 
1 pve a® ; ; 
Then there is the effect of the additional terms (14) which, however, only 
occurs for terms of O(R-5) in the series (6). 

The asymptotic series have now been found for the case of unsteady 
flow in which no simple exact solutions can be found. A feature of this 
asymptotic solution is that conditions at time ¢ depend only upon the 
value of R and of v, and its time-derivatives at that time and not upon 

U) previous happenings. This corresponds to the result of steady boundary- 
layer theory that conditions at a point where the suction velocity is large 
depend only upon the values at that point of the suction velocity, the 
stream velocity and their derivatives along the surface. The present 
results therefore not only will provide a valuable check on coming experi- 
mental work but also provide interesting comparisons with conventional 
boundary-layer theory. 

) 
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THE EQUILIBRIUM DISTRIBUTION OF 
THE LONG-PERIOD TIDES OVER AN OCEAN 
COVERING THE NORTHERN HEMISPHERE 


By M. R. REES (The University, Reading) 


[Received 12 July 1949] 


SUMMARY 

tesults of a calculation of the equilibrium distribution of the long-period tides 
over an ocean covering the northern hemisphere of a uniform spherical globe ar 
given here making allowance for the change in the gravitational field produced by 
the tides themselves. The calculation is based on a normalization process described 
by M. Brillouin in 1927. The main result shows that, on the assumptions on which 
Brillouin’s paper is based, making allowance for the mutual attraction of the elevated 
water gives an elevation which approximates very closely to that which would be 
obtained by multiplying the uncorrected level by the constant 1-126. 


1. Introduction 

By the equilibrium distribution of tides is meant the distribution as 
calculated from the tide-generating forces of the moon or sun on statical 
principles. The problem of making such a calculation over an ocean 
bounded by continents is simple so long as the mutual attraction of the 
elevated water is neglected. The solution is given in Lamb’s Hydro- 
dynamics (1). If there are no continents, so that the ocean covers the 
whole globe, then the problem of allowing for the mutual attraction of the 
elevated water is also simple, and the solution is given by Thomson and 
Tait (2). But when both continents and elevated water are taken into 
account, the problem is difficult, and no case has previously been worked 
out in detail. A general analysis of the problem was given by H. Poincaré 
(3) in 1895, and another by M. Brillouin (4) in 1927. The present paper 
consists of an application of Brillouin’s method to the long-period tides 
in an ocean covering the northern hemisphere. 

I wish to thank Professor Proudman, F.R.S., very sincerely, for suggest- 
ing this problem and for giving me much valuable assistance. I am also 
indebted to Dr. P. White for pointing out a mistake. 

2. Notation and introductory theory 
(r,@,), spherical polar coordinates with origin at the centre of the 
earth which is assumed spherical; 7 denoting the radius, @ the co- 
latitude, and ¢ the longitude. 


a, the mean radius of the earth. 
y, the constant of gravitation. 


[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 
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where ds (= asin @ déd¢) is an element of the spherical surface and all 
the integrals are taken over the whole sphere. 


The value of the constant D in V, is calculated by using the fact that 
Cds 0. 
where the integration extends over the ocean. 


3. Long-period tides in an ocean covering the northern hemisphere 


For an ocean bounded by the equator, 


k (4rrypa)/g 3p/p if O<O0< 72/2, 
0 if /2<O0<7. 
The effective potential producing the long-period tides does not depend 
on ¢, so the coefficients A, of § 2 are zero except where n is a perfect 
square, that is when S,, is a zonal harmonic of order vn. Consequently, 
S,, may be replaced by P,(cos@) and ¢ may be replaced by » everywhere. 
- It follows that 


: ] : ; . . _ 
t C [Con Fo Cyn Tat. +Cy-in 74 | (kk—2n 1)P..|, 
and on substituting for ‘YY (r i n—1) and rearranging 
_— 
A C, [Bon(h—1) Py +64, (k—-3)P, +... +0, n(k—2n+ DP, | 
Mo,,(k 1) Py a, ,(k 3)P, ae a, ,(k 2n- 1)P, 


where the a,,, are constants and C,,,,, C,.,, are now calculated from 


nm? myn 
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™ Since | P,P, sin@ dé = 0 if s and n are both even or both odd, 
hat é ™ 
l Me] 3 S— 9), Bias ae. « , 
*:***"" if s is even and n is odd, 
S | ($/2 1) 2}! 
.d | PP. siné dt | PP. siné dé, 
here ) 7 
t follows that 
i 
Qa (2n P| S a Zé 1)P. sin 6 dé (mm n) 
Ls=90 
pel 
rfe —~s a (ns 1) —] 1 )( ,) 2] .3....(8 1).1 3. n 
nth ay s—n)(s+n-+-1) (s 2)! |(n 1)/2]! (2) 
lere, if s is even and v7 is odd 
. ” 112 l j\s+n—-i2] .3....8.1.3....(n—1) 
 —* . s)(s+n- 1) |(s 1)/2|! (n 2)! 


if s is odd and v is even. 
The calculation consists of three parts which are here labelled the first, 


second, and third process respectively. In the first process, using the 


juations (1) and (2), the coefficients C,, ,, and C,,,,, and the functions YF, 
re calculated as far as 8. They are calculated in the following order: 
w. 
0,0 Wo 0, Fo; Con (0 - n) 
: , 
C4; 0; © nH; GQ, (l<n) 
. o yw . Y 
Lit: fi: Tt Ga Oe 


id the results are tabulated 

In the second process the coefficients A, (n = 0, 1,..., 8) are calculated 
nd tabulated, and the expression for the tidal elevation given in terms 
f Legendre polynomials and the undetermined constant D. 

In the third process the constant D is evaluated by using the fact that 
e total amount of water remains constant, and the tidal elevation is 
ulculated at co-latitudes 0°. 10°..... 90°. For an ocean covering the whole 
globe the allowance for the mutual attraction of the elevated water can 
made by multiplying the uncorrected level by the constant multiple 
1-126(1). The tidal elevation is now compared withthe uncorrected level and 
vith the uncorrected level multiplied by 1-126 at co-latitudes 0°, 10°,..., 90° 


nd the results are tabulated 
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The results of the calculation show that, for an ocean covering the 
northern hemisphere, a good approximation to the actual height is given 
by multiplying the uncorrected level by the constant 1-126. 


k irst proce SS 


| 
| 
For convenience. the calculated expressions for Q.; ®M.. and _ a for | 
n | ao 8, are collected together in Table 10. The appropriate values 
for a, ,, in equation (2) can therefore be read off from this table remember- | 
ing that 4, , = Pgin/Cn» Substituting those values in equation (2) gives 
Cin.n Lor m n. Equation (1) gives C’, , for a eee 8. The results are 


given in Tables 1 to 9. 


TABLE ] 


TABLE 2 


TABLE 3 
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the TABLE 4 


iven 


C3.n/C3.3 
I 
2 0489 
7 Io14 
1359 
lOY 
O*O1 sol 
lues S 1374 
bey 
ive 
TABLE 5 
i 
, tn) Can/Cas 
I 
$i 430° 
{ 06330 
if oo1055 
1259 o*o0C $95 
TABLE 6 
( ( 
I 
2115 
7264 7476 
34 
[ABLE 7 
Cen/Ce 
I 
54 
15054 Cc 21532 
TABLE 8 
( C e¢ 
I 








Values of ¢ 













Onr=a, C,,90 ( 
Second process 
| 
1338 
3254 
T23 
Values of A, 


go 
gH (P,+-D)+-gl 
-[0-04285 
-|0-03802-4 

+ 


+[0-006440+ D( 


0-01163 


Tih. ;' 
1 a Ends, Va 


gH (P,4+-D)+ ¥ 


D( 





,enad 0, T (a= Q I... 8): 
n n = 
SO = 56,8. ¢,, Fix Disk Lipe. 
r 0 
TABLE 1] 
] Ht 
) 1 ) 
) I ay ) 
7) ay 1214) 
-\) ~ By 5020) 
1) Di 3450) 
) Di ) 
) 1 Dy 1569) 
) | Dy 1s67) 
) 2 Di 5704) 
gH and A, /qHC,,,, (n a 8). 


A 


10) 


n n 


(using Table 12) 

/{|0-01218-++-D(0-4150) |P, 

D(0-2170) |P, +-[0-06180-+- D(0-006668) |P, 
0-05102) |P;+-[0-00039-+ D(—0-002119)]P, 

+ D(0-02516) |P;+-[ —0-00025+- D(0-001109)]P, 


0-01559) |P,-+-[0-000242-+ D(—0-000870) |P,}. 
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Third process 
Since ( ( € ds = 0, where the integration is taken over the ocean, THE 
D —0-01819. TO va 
C/H = P,—0-01819+-(0-01218—0-00755)P, , 
+ (0-04285—0-00395)P, + (0-06180—0-00012)P, 
+ (0°03802-+- 0-00093) P; + (0-00039-+- 0-00004)P, 
+ (—0-01163 —0-00046)P, + (—0-00025—0-00002)P, 


+ (0:00644-+- 0-00028)P, + (0-00024-+-0-00002)P,. 


TABLE 13 











a onan , The p 
(] ° ro° | 20° | 30° | 40° | 50 panied | 
: Re i Zz ‘ae that this 
CH ['121 1:069 | 0920 | 0°696 | 0°423 | 0132 | 
4) H I yee 0824 of 0280 o'120 ind the 
€,./H | 11 1075 | 0-928 "704 | 0°428 | 0135 high-pre 
—— — —- - ——_—_—— As usua 
arn penne that ra 
6 ¢ 70 &¢ 90 “8 wk 
= ; ise . uniform 
: H or144 | 0° 266 0504 | 0°544 wave m 
G/H | —or125 | —o-324 | —o-455 | —0'5 With 
¢,/H O141 | —0°365 51 0°563 
Ae, a ee | MAS ae aon Te wise worked 
Table giving the tidal elevation at co-latitudes 0°, 10°, ..., 90°. flow wi 
approx 
In these tables f denotes the tidal elevation when allowance is made for compat 
the attraction of the elevated water, ¢, denotes the uncorrected level, owiel. 
‘ee cae to be n 
¢, denotes the uncorrected level multiplied by 1-126 
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THE APPLICATION OF BOUNDARY-LAYER THEORY 
TO SWIRLING LIQUID FLOW THROUGH A NOZZLE 


By A. M. BINNIE and D. P. HARRIS 


(Engineering Laboratory, Cambridge) 
[Received 27 September 1949] 


SUMMARY 


The passage of a swirling liquid through a convergent-divergent nozzle is accom- 
ed by the formation of a retarded layer on the nozzle wall. On the assumption 
this layer is thin, the motion in the main stream may be taken as frictionless ; 


the pressure and velocity distributions in it, when the flow is derived from a 


pressure reservoir, can be obtained by an application of critical flow theory. 
\s usual, it is necessary to suppose that the angle of the conical nozzle is small so 
radial velocities can be disregarded. The streaming velocity is found to be 
niform over each cross-section, and it is equal at the throat to the velocity of a ‘long’ 
moving on the surface of the air core. 
With the aid of the results mentioned above, the boundary-layer equations are 
rked out by an extension of the method used by Taylor in his study of swirling 
w without streaming, and a solution is obtained with the use of the Pohlhausen 
proximation. A numerical case is examined in detail, and the boundary layer is 
for mpared with those caused by swirl with no streaming and by streaming with no 
el. virl. The effects of surface-tension over the air core are considered and are shown 
negligible in the numerical example. 


1. Introduction 
(HE problem of swirling fluid motion inside a conical surface was attacked 
by Taylor (1). He considered a state in which the axial velocity in the 
iin stream is negligible compared with the swirl; and, on applying the 
\pproximations of boundary-layer theory, he found that the retarded layer 
ear the wall cannot maintain its position but is forced towards the apex 
f the cone. The results were pushed to a numerical conclusion, and in 
ypical examples the boundary layer was shown to occupy a very large 
lraction of the cross-section at the outlet end of the conical frustum. 
Thus in an actual nozzle it may have an important influence upon the 
ischarge. Tay lor’s work also explains the action of the “¢ tyclone’ method 
{ dust extraction from air, which was tested by Linden (2). In this 
levice the contaminated air is admitted tangentially to a conical vessel. 
The particles of dust are impelled by centrifugal force to the curved wall, 
long which they move in the boundary layer towards the closed apex end, 
vhere they are trapped. The cleaned air escapes on the axis at the other 


end of the vessel 





[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 
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The swirling flow of a perfect liquid through a convergent-divergent 
nozzle was examined by Binnie and Hookings (3) and by Binnie (4), who 
obtained an expression for the discharge under given conditions. Making 
Reynolds’s assumption that the radial velocities can be neglected, they 
proved that the axial velocity is uniform over the cross-section; with 
unimpeded flow, it attains at the throat a value equal to the velocity of 
a ‘long’ wave moving axially on the spinning surface of the core. Their 
analysis deals with motion under gravity, but it may easily be modified 
to apply to the case, now to be considered, of flow derived from a high- 
pressure reservoir. 

In the following pages an attempt is made, in the manner described by 
Taylor, to derive the laminar boundary-layer equations when the main 
stream of liquid possesses both swirl and axial velocity, the latter being 
calculated by means of the approximate method mentioned in the previous 
paragraph. The equations are then solved in order to provide estimates 
of the thickness of the layer and of the discharge through it. 


2. The motion in the main stream 
We consider a nozzle discharging freely and supplied from a tank in 
which the absolute pressure is P, measured at a point where the velocity 
is negligible. The flow being irrotational, the swirl velocity V at radius r 
is given by an expression of the form 
V = Qryr, (2.1) 
in which © is constant throughout the liquid. At a cross-section where 
the radii of the nozzle and of the core are a and 5, it follows from Bernoulli’s 
equation that if W and p are the streaming velocity and the absolute 
pressure in a liquid of density p, 


Q2 ) 
P | 5(2+ 5) f , (2.2 
p 2 r2 p 


The radial velocity here has been omitted as being very small compared 
with the other two components. Now the variation of p over the cross- 
section is seen to be 


Pe 3) LY (2.3) 
p 2\62 #2 
hence from (2.2) it appears that W has the uniform value 
9P QA! 
W = (I. (2.4) 
po 
The discharge Q is then given by 
Q = n(a?—b?)W. (2.9) 
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gent | The axial coordinate in the nozzle is denoted by z, measured positively in 


who lithe direction of streaming; and on differentiating (2.4) and (2.5) with 
king | spect to z, and eliminating db/dz, we find that 
they . ’ - 5 
Q?a da Ww?—w? du 
“47 1 (2 6 
W1Th i ~ ie ) 
h dz \ dz 
J ol 
’ y2(772 2 1 
heir | ner ‘« (2.7) 
A | Dh 4 | 
lhed ys 
igh s shown by Binnie and Hookings (3) to be the axial velocity of a ‘long’ 
ve on the free surface. Now at the throat, da/dz = 0 and dW/dz + 0, 
L by erefore at that cross-section VJ WW, and from (2.4) and (2.7) a quadratic 
lain uation for 6? is obtained which yields 
i ' _—— 
5 bf pQ* | 16 Pa7z\} 
ous ; L+|] iy) (2.8) 
a; 8 Pa? \ p&2- | 
ites 
suffix ¢ referring to the throat. The other root, being negative, is 
relevant. From (2.5) and (2.7) we deduce, as the final expression for Q, 
QO r 2/az 
| a (2.9) 
: mn Qa, V2 b?/a? 
| in this (2.8) is to be substituted. The non-dimensional quantities 
IS i ; 
olved are shown in Fig. 1, where Q/(Qa,) and b,/a,; are plotted on a base 
I o92/(8Pa?). The discharge becomes zero when pQ?/(8Pa?) = }, i.e. 
?.] hen $0?/a? P/p, and it vanishes because the cross-section of the stream 
ere lisappears and the whole initial pressure is turned into swirl, leaving 
li’s thing to produce streaming velocity. At the other extreme when Q — 0, 
ate tis found from (2.8) that b?/a? Q/a,)/(P/p), and from (2.9) Q becomes 
702(2P/p)' as it should. Fig. 1 also shows that the ratio b,/a,, if small, is 
ery sensitive to slight alterations in 2 and P. 
2) For insertion into the boundary-layer equations we shall want the 
ilues along the nozzle of W, dp,/dz, and dW /dz, p, being the pressure at 
- the wall. When Q has been determined, W and 6 may be removed from 
SS- 2.2), (2.3), and (2.5), leaving an expression for p, in the form of a cubic 
equation; and after this is solved, W is calculated from (2.2). An easier 
9 . 
” procedure, however, is to eliminate p, and b from these three equations, 
r the resulting cubic equation for W, viz. 
Q i, Fias P Q 
W3 YW? | 2—|\W+2 ¥ = 0, (2.10) 
+) 7a" L“ p p 7a* 


is simpler coefticients. One of the roots of (2.10) is negative and therefore 
irrelevant. The other two are positive, the smaller being the velocity 


upstream from the throat and the larger being the velocity below the 
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it at the section having the same area; at the throat these two roots 
equal. The velocity gradient is obtained by the differentiation of 


10) which yields 


2 da Q ( .. PO 
— We W +-2— A 





| J > ” 
1} a dz T¢ g* p Tra 
(2.11) 
2Q) 0)? P 
ai — Vj 2 
ra a* p 
throat this expression assumes a 0/0 form, because there da/dz and 


denominator are zero. On applying the usual method for dealing with 


LS we find th the aid of (2.7) and (2.9) that 
di Q( da/dz2 \ (2.19) 
dz a,(2-+-b?/a?)} _— 


s the radius of curvature of the wall, which is the reciprocal of d?a/dz, 
st be known at the throat in order that the velocity gradient may be 
termined. Finally the pressure gradients are calculated by differen- 


» 9° 


he theory explained above is evidently subject to certain limitations. 


when the angle of the cone is large the assumption of negligible 
lial velocities is violated. Secondly, it is a matter of common observa- 
that in actual nozzles the discharge at high swirls does not cease: 


ler these conditions the whole of the discharge must pass through the 
ndary layer. Measurements made to elucidate this point were shown 
sinnie and Hookings (3) in their Table 1, from which it appears that 
simple theory does not break down if the swirl be moderate. 


\s a numerical example we have considered the nozzle indicated in 


rig. 2, The convergent part consists principally of a cone of inlet radius 
2-59 in. and semi-angle 10°, extending for an axial length of 10-87 in. To 
luce parallel streaming motion at the throat, the nozzle wall is con- 
ied by a succession of conical portions having semi-angles 8°, 6°, 4 
| 2°, each of axial length $1 lo form the throat itself, a circular are 
placed tangentially to the 2° cone; equation (2.12) requires that the 


lius of this are should be specified, and it has been taken as 1 in. The 
radius a, is } in., and the nozzle is symmetrical about the throat. 
[he nozzle operates with water under a head of 30 ft., so that 


Pp 30 x 32-2 ft.?/sec.2, 
{ irrotati constant £2 0-7 ft.2/seec. Hence 
oQ2?/(8 Pa?) 0-0365, 


d it is seen from Fig. 1 that the motion is far removed from the range 


ear pQ?/(8 Pa? where the theory is known to be unsatisfactory. 
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The calculations have been pursued for a short distance downstream | the ap] 
from the throat, and the values of W, dW /dz, and dp,/dz, which are required | jrregulé 
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Fic. 4. Distribution of streaming velocity and wall pressure gradients in 
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for the solution of the boundary-layer equations, are displayed in Figs. 3 
and 4; in the former the values of p, are also inserted, and the shape ot 
whe 


the air core has been added in broken lines to Fig. 2, which shows that in 
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tream| the approach to the throat the core diameter is almost constant. Slight 
puire rf cularities are visible, due to the absence of smooth transitions between 

e various conical portions of the wall; but they are not important, since 
e theory is approximate and, moreover, the boundary-layer equations 
n only be solved by a step-by-step method. The computations were 
ind to involve small differences between large quantities, consequently 
euse of a calculating machine was necessary to secure sufficient accuracy. 


Near the throat the results are evidently sensitive to the form of the nozzle 
ere. hi 


pect the flow to be modified, for the theory set out in this section does 


a nozzle terminating abruptly at the geometrical throat, we must 


contemplate a rapid reduction to zero of the wall pressure in this 


ion 


3, Effects of surface-tension 


We must now consider whether the forces due to surface-tension have 


ppreciable effects. In section 2 the pressure over the surface of the core 
s taken as zero, but when surface-tension (denoted by y) is taken into 
unt, this pressure is 
p y( 4 ) (3.1) 
RK, R, 


ere R, and R, are the principal radii of curvature. Since the free 
face is a surface of revolution, R, is the intercept on the normal 
tween the surface and the axis, and R, is the radius of the generating 
ve. We will confine our attention to the neighbourhood of the throat, 


that R, may be taken as equal to b and R, as very large. Accordingly 


z 3 ] reduces tO _ b. (3.2) 
aid ) Q)-/ ] l y 

nd (2.3) beeomes l : ) Rony (3.3) 

p 2 \b? pb 
ernoulli’s equation (2.2) is unchanged, and when used with (3.3) it yields 

2P QF 2by 
We {i cay (3.4) 

0 b? pd" 
ving that the streaming velocity is again uniform over the annular 


ss-section. Employing the same procedure as before, we differentiate 


2.5) and (3.4) with respect to z, and, after eliminating db/dz, we find that 


= : (3.5) 


()-a by \da W? Wedw 
| pQ | lz W dz 


(3.6) 
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The flow under consideration is seen to be impossible if by + pQ?; and, 
subject to this restriction, it follows from (3.5) that at the throat IV W. 


re 


[t may readily be proved that W, is equal to the velocity V’ of a ‘long’ 
annular wave moving axially on the surface of the spinning liquid, 
Adapting the method employed by Binnie and Hookings, we take q to be 
the additional velocity due to the wave motion at points where the 
infinitely small surface elevation is n. The equation of continuity is then 


‘7 a2—bh?) 2aby\(V’ q) (a= b?)V’. (: 


> 
v0, 


With surface-tension operative, the increase of pressure at any radius due 


to the wave is 


se l l Y| l ) : 
‘ ; (3.8 
2 (h n)* ‘a p h UW] bh 
: (2? by 
which reduces to " ( A (3.9 
b3 p&2- 


when 7 is small. Hence, dp being the excess of pressure due to the wave 


Oo ()*, Dy - > v9 
p , & (1. ”) 1(V’ +9)? = 1V"2, (3.10) 
p b3 \ p&2* 7 x 


/ 


motion, 


On substitution from (3.7), (3.10) becomes 


op 0)? n l- by } 2byV"? (3.11 
p b3 pQ2d*, a2—h2’ 


and the surface condition 6p = 0 is satisfied when V’ has the value 
assigned to W, in (3.6). It therefore appears that the general conclusions 
of critical flow theory are not vitiated by the action of surface-tension. 

We next attempt to determine b, by equating the throat values of (3.4) 
and (3.6), but in place of a quadratic equation for b? we now arrive at the 
biquadratic equation 


l =e , a ° 
b} 3A b? b? Ka? b, az 0. (3.12) 
2M 
p62" , Y ‘ 9 
where M and A fas (3.13) 
8P pid” 


This may be solved by trial, a general solution not being obtainable 
without serious complications. Some immediate progress, however, can 
be made if A is assumed small, so that the solution does not differ widely 
from (2.8). When K is zero, (3.12) reduces to the form previously used, 
of which the root, say b, = 8, is given by (2.8) and shown in Fig. 1. We 


therefore suppose that — (3.14) 
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und here € is small, and, on inserting this into (3.12) we obtain approximately 
We. 


3p°-+-a (3.15 
ne | ? (28/K)§(B2/M)—1}" — 


“ro 


'| For water, the ratio y/p (which may be termed the kinematic surface- 
tension) is about 0:00261 ft.3/sec.2. In the numerical case under considera- 








ion it is found from (2.8) that 3 = 0-02317 ft. The corresponding 
hel 1 . 9 @e 
ue of By/(pQ*) is only 1-24 10-4, and we see from (3.4) and (3.6) that 
3.7 e effects of surface-tension are negligible. This conclusion is confirmed 
’ 3.15) which yields « 0-52 x 10-4. 
4, The motion in the boundary layer 
i It is now necessary to employ spherical polar coordinates. The notation 
sed is set out in Fig. 5, which indicates how the component velocities 
oe er 2 i of 
ees -— ae = 
a er ieee 
‘f- 
1] 
vlue 
Fic. 5. Svster fs} rical polar coordinates. 
lO 
, w are related to the axes. For flow symmetrical about the axis 0 = 0, 
3.4 the complete equations of ste LcLy motion for a viscous liquid are proved 
the Goldstein (5) to be 
Cu vecu Uu pe l op — 2u 2 av 2vcoté (4.1) 
A vi V-u te heen mal a ¥ 
( R R oo R R p R | R R? Cc g k- ‘ 
Cow » OW wu vw cot " Ww 
u = vj V?w—-——— }, (4.2) 
oR’ Reo R R R* sin?6 
13 : . : : " 
ere p is the pressure, p and v are the density and kinematic viscosity, 
hl bi l ¢ ] C ‘ C ? 
— 72 [ R aa sin @—),. (4.3) 
cal R? R C R R? sin 6 06 fale 
lel) The complete equation of continuity is 
sed 9 : 
Cu ou lL co v 
We ,coté 0. (4.4) 
ci ki R Rc6 R 
l4 


Many of the terms in these equations are small, and, as shown by Taylor, 


H 
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the forms of (4.1), (4.2), and (4.4) which may be used in an approximate | becaus' 


solution are 








cu vou w 1 dp v ou , 
U “i + a net tok eo (4.5) 
oR R 8 R por R? 0G? 
ow. vow. wu v Ow (4 
o— + + “ore ee .6) | When | 
eR’ RO’ R RCC )y 
mome! 
ou, gu, 1 & e 
-+2—+ Q. (4.7) 
OR RR’ Roe ‘ 
We assume that the boundary layer is thin and that the pressure “s 
’ u—<¢ 
through its thickness is the same as in the adjoining main stream. It was 
shown in section 2 how an expression for this pressure could be derived, 
but it is too complicated to be of service here, and for the time being we 
shall take the pressure to be given by : 
: The 
lop : =? 
Pa F(R), (4.8) | 2 relat 
p¢ R : 
where F(R) is the function of R calculated from (2.10), (2.11), and (2.2). 
For the same reason we suppose that uw in the main stream close to the 5 
boundary layer is of the form 
u G(R). : (4.9) 

The boundary layer is taken to have a definite thickness 6 which is a 
function of R, and the momentum integrals are obtained by integrating | !0r < 
the equations of motion (4.5) and (4.6) between the limits 6 v at the | andv 
wall and @ = «—8/R at the edge of the main stream. In this way (4.5) and | 4 

{ 
(4.8) yield 
: : : “8R 
cu vcUu . we vicu 
nu do d6-- | k do (4,10) 
; oR J Ret J \ R R?| c6|, which 
r—0/R x—0/FR x— 0/R We 
eu/e6 at 0 .—6/R being supposed zero. From this equation v can be 
removed with the aid of (4.7), for 
XY Y | 
- wv eu 1 (e(uv) ov so tha 
do= | uw" \ dO : 
R e0 R\ © aad ound 
r—0/R xr—d/R meth« 
x 
l ¥ cu . 2u? sisten 
[wr set | (uaa +—) a0. (4.11) 
R 1-3/R oR' R in the 
r—0/R 
To evaluate v,_5;, we integrate (4.7) which yields 
“s and 
cu 
Va-S/R | : 2n) dé (4.12) 
) i R 
= where 
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ecause v, 0, and thus we find on using (4.9) that 
‘i ou ‘ 
uv G (R + 2u) dé. (4.13) 
la—O/R \ OR - 
O'R ; 


When (4.13) and (4.11) are inserted in (4.10), the result is the longitudinal 


omentum integral 





f- # a ary 
7 ~U : dé 2 : dé } R‘ > ~ n= 2u dé i 
p c rf 2 t t . \ Cc R 
] O/ ER i 
"ofp, sa i LIS 4.14 
‘y's R|20),° 
r—O/R 


The same method can be applied to (4.6) which, with the assistance of 


relation analogous to (4.11), yields 


cw l r os l 2uw 
Ti ' dé pul | (are 4 dé + 
( R R : fe J cR R 
R 
~ wu v |ow z 
dé : , (4.15) 
R R?| 00 |, 
. 0=a 
x OR 
6 at @ .—8/R is taken to be zero. At this limit w Q/(R sin «) 


dvis given by 1.12), henee t.15) becomes 





" (uw) Q) ¥ u ww v lew 
dé (R 2u} d0 3 dé = 
F ( R R SI a J ( R J R R?2 o6 r) 
i / r—o/R 
(4.16) 
hich is the transverse momentum integral. 
We now use instead of 6 asa variable, ” being defined by 
R 
(ax 0) (4.17) 
that it varies from 0 to 1 as we pass from the wall to the edge of the 
indar\ layer \t this Stag Is necessary to adopt the Pohlhausen 
thod of approximation, in which simple arbitrary assumptions, con- 
t with the boundary conditions, are made for the velocity distribution 
the bound ur'y lavel Accord oly we suppose that 
G R)d(n) (4.18) 
Q) 
(7), (4.19) 
ki SIL a 


(4.20) 
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It will be seen that (4.18) and (4.19) satisfy the boundary conditions 
u Ww 0 at 7 O; u G(R), w= Q/(Rsina), eu/éen ew/on = 0 
at 7 A 

Before the integrals (4.14) and (4.16) can be evaluated, @u/@R and 
O(uw)/OR, which are differentials at constant #, must be expressed in 
terms of ¢, 6, and 7. We have from (4.18) that 


(= rp ps ll (4.21) | 
C R A dy c R fa) dR 
but from (4.17) 
on .—G uy, 8) ds n 1 R déd\ (4.2 
- — (x — — : 22) 
C R A a) rote dR R ra) aR} j 
4 dd Y R dd ' dG 
a (—*) G— 1] Lt. 4,23) 
ae @R)y ‘dy R\ 8 dR)'*aR ae 


Similarly, from (4.18) and (4.19) 


sin Y = aed _2G¢ dd 21 R dd 1 2 l dG G ; (4.24) | 
Q oR }, R dy k 5dR RdR_ R 


We now change the variable in (4.14) and (4.16) from dé to —(8/R) dy 
and utilize (4.23) and (4.24). A simplification is effected by means of the 


expressions 1 


1 
- ] 7 19 
nd- ¢ dy 1 $ dy (4.25) | 
dy 7 I 2 
0 0 
; ] : 
and n- p dy ] | d dy, (4.26) § 
dn ‘ 
0 0 


which are obtained on integrating by parts and making use of the boundary 


conditions. The momentum integrals then reduce to 


a 1 1 
G  d@ G déd*\/ : dG Q?2 , F 
Bel Mai | ddan | $ an) a | $? dy —— 
Ri dk 28°dR}\, 4 dR GR sin*a) . G 
0 0 


0 


an le . (4.27) 
0" dy. n=0 


1 1 
G dG G ds? i. ¥ Vv dd 
l L 1 | gd Oe hn 4.28 
(i "dR ae | ile be in) mal. : oe | 
0 0 


These may be put in non-dimensional form by substituting 


t PS Q 1 
—- 5, | — Jt (4.29) 


R,’ Ry\v sin « 
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and replacing (4.8) and (4.9) by 


Pa _ F(R he 
p¢ R — Ro sin2a! ' ) 
Q) 
nd u G(R) g(R,), 


Ry sin a’ 
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(4.30) 


(4.31) 


vhere R, is the distance along the’ generator from the apex to the base 


fthe cone. Then (4.27) and (4.28) become 


l 
(9,4 | § 
\R 


0 


0 


1 
J aos \I "49 dg ats f 
; @ dy d- d , EO ; | 2 1 oe 
dR, 20, dk, | | ] 7 dk, al 4 ea J 
U 0 


] Fa 
oF dy n * 


1 1 
qd dq g dd; r ; — L{dd 
3 bdr | ‘deh =i | 
| R, dR, 205 aR,)| | , p ii ) oF dy n=0 
0 


(4.32) 


(4.33) 


\part from boundary conditions we have not yet employed the assumed 


expression for ¢ indicated in (4.20). With this particular form 


1 1 
: , ' ‘ d 
p* dy 0-53. db dy 0-6, | ‘| 2; 
A dy n=0 
0 0 
in that case (4.32) and (4.33) reduce to 
qd dd; 6 dg 9 g 5 lof 30 


) 2 + -_— -- 
‘2 AR 3 52 ? 
dj dh, LR, Rk, gk} g oF 
gd dd? 4 dg “ g ; 15 
252 2 : —* 
26¢dR, dR, R, 3: 
If we replace (4.20) by the more elaborate expression 
d 2n 2n°- 7}, 
vhich also satisfies the boundary conditions, we find that 
1 l 


d2 di 0O-5825. db dy Ord, 
. dy. 7=0 
0 0 


and that (4.35) and (4.36) are changed to 


( dd- dé 9-9 5 j . - 
bot an FREE 9 199 , 9915 _ i7.99F — 
OF dk, dk, R, gh} q OT 
g di. dg . g 17-02 


22dR,'dR,'R, 8&8 


: bs ; 


(4.34) 


(4.35) 


(4.36) 


(4.39) 


(4.40) 
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It will be seen that the alterations in the coefficients are not large. Some | 


idea of the magnitude of the approximations can be obtained from the 
subtraction of (4.27) from (4.28) which yields 


aG@_ oD pli gag 
7 -—_—_—- +f / | ¢*dyn = 0. (4.41) 
dR RF sin?« /I l 
0 
The coefficients of F in the two cases considered are 1-875 and 1-717. 


For comparison an expression, which is exact as far as this section of the 
work is concerned, can be derived by differentiating (2.2) after putting 
It then follows that 


r a. 


dG Q? 


G = s 
dR Rsin?a« 


F = 0, (4.42) 
indicating a discrepancy in the final term. 

Evidently the solution obtained above is not valid when the streaming 
velocity is zero. In those circumstances (4.18) shows that « would also be 
zero, whereas Taylor's theory yields a definite expression for this velocity. 
However, at the other extreme, where Q — 0, a check on the work can be 
obtained by comparing the results with those due to Mangler (6) and in 
the next section the metliods used above will be applied to this case. 


5. Boundary-layer theory when the swirl is zero 

We are now concerned with simple flow without swirl through the 
nozzle. As a first step it is necessary to determine the conditions in the 
main stream, using the approximations and notation of section 2. Taking 


the throat pressure to be zero, we have 


)»P\l 
Q nai J’ (5.1) 
Pp 
2/9 P\1 
and it follows that i} a= }" (5.2) 
a- P ] 
> 1 
whence from (2.2), Pa i — : (5.3) 
pp at 
Accordingly (4.8) and (4.9) are replaced by 
l ) 5 4 > 2 
Pe ac = (5.4) 
pc R 2° sinta p 
2 9P\1 
and u G(R) a“ & )* (5.5) 


R? sin2a 
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Equation (4.6 disappears; and as (4.5) leads to (4.27), in the latter we 
substitute (5.4) and (5.5) and put Q 0. The result is 


1 1 1 
R dé ; ! oer v Rsin2a Id 
= | | ddn— | ¢?dy)—2 | 2 dn+2= 54 1. 
25° dR ; . 8° a7(2P/p)*| dy |,-0 
0 0 U0 
(5.6) 
lo get this into non-dimensional form we write 
R a?(2P/p)*\4 ee 
R, : Os 8) He . p) \2. (5.7) 
R, ; | v.R3 sina | 


then (5.6) reduces to 


Os * >: 
= iB | | g dy | p° dy 2 | d” dy 19 Ry 


0 0 


ad os 
. (6.8) 
dy n=0 


With the numerical values given in (4.34) this equation for 63 becomes 


> 


05 





dds ey ‘ - 
TR. 16 R 30 R35, (5.9) 
sted: a 


e solution of which, subject to the condition that 5, Oat #, = 1, is 
82 0 ¢3(] Ri). (5.10) 


The term R} is appreciable only in the initial stage, and over the rest of 


the cone the momentum thickness is 
uv (35) Ry } d db") dy 0-203 R}. (5.11) 


Mr. E. J. Watson has pointed out to us that, for comparison with 
5.10) and (5.11), use may be made of Mangler’s work (6) on exact solutions 
f axially symmetrical boundary-layer equations. As a special case of his 
general theory, Mangler considers radial flow over an unlimited flat plate 
towards a central hole, and since the velocity in the main stream varies 
nverseiy as the square of the distance from the origin, his analysis applies 
lso to a cone of small angle. His boundary conditions are different from 
ours, for both at the origin and also at an infinite distance from the wall, 
he takes the whole of the boundary layer to have the same tangential 
elocity as the main stream. It can then be shown that the momentum 
thickness is 0:291R;. The agreement is closer when displacement thick- 
nesses are compared; (5.10) then leads to 0-506 Rj, while Mangler’s result 


IS 0-587 R;. 
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6. Numerical example of boundary-layer calculation 
For the purpose of step-by-step calculation (4.35) may be written in 
the form 























—; so. oe 
—_—— + A OF . (6.1) 
dk, g 
if ( Ss i 
where A= i 6 “og +72 Sil A — is! (6.2) 
g dR, Rk, gk} g 
[ | 
(q) - aii i 
o:O6 r— - a erie LL ee oe a - _| 
\ 
\ 
~\ 
| | - +. 
| | (b) . Vt 
ail \t 
N | L 7 \! 
+ — —i 
| 
TANCE rROomM Ss THROAT iN 1O N =HES | 
T T 
oO-+ os o Ss Oo 
R= FYR, FOR 10° CONE 


Fic. 6. Thickness of boundary layer. (a) With swirl only (Taylor’s result). 
(6) With streaming only (equation 5.10). (c) With both streaming and swirl 
(equation 4.35). 
is assumed constant over each step. The coordinate z of section 2 is 
identified with — R of section 4, and to evaluate g and A use can be made 
of the results displayed in Figs. 3 and 4. Denoting the ends of consecutive 
steps by the suffixes m and n-+-1, we then find on integrating (6.1) that 


(82) ,,43—30/(9n 4A) - 
~ = exp] —A{(R,),,.,—(R,),} 1. (6.3) 
(5?),, —30/(g,, A) pl 1 the 1 ( wnt 
The process was carried out starting at the point 6, = 0, R, = 1, adjust- 


ment of the origin of R, being required at the sections where the cone 
angle alters. The outcome is shown in dimensional units by curve (¢) 
in Fig. 6, in which the ends of the steps are indicated by dots. At the 
throat the boundary layer has a thickness of 0-00695 in., and it occupies 
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, fraction 1/24-85 of the cross-section of the stream. For purposes of 
en i mparison two other curves, extending over the 10° portion of the nozzle, 
ive been added to Fig. 6. The upper is Taylor’s result for zero streaming 
(6.] elocity and the same swirl Q 0-7 ft.2/sec. as in the motion now under 
consideration, and the lower is calculated for zero swirl from (5.10), the 
reservoir pressure being chosen to produce the same streaming velocity at 


the throat as in our case. 
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" Fic. 7. Reduction of discharge due to boundary layer. 
In most boundary layer calculations which have hitherto been carried 
ay ° . . ° ° — ° 
av out, the field occupied by the fluid has been partially infinite in extent. 
Here, however, the field is enclosed, and the condition exists that the total 
e discharges over all cross-sections of the nozzle must be identical. It is 
2 is a ; er hap he 
: therefore desirable that the reduction in discharge due to the retarded 
nade 1 ‘ ° vs . . 
' layer should be estimated and compared with the frictionless discharge Q. 
utive wn 
This reduction is given by 
hat 
1 
(6.3) W2zad | (1—d) dy aL 70, (6.4) 
just- ind is plotted in Fig. 7 as a percentage of Q. It is seen to be very small 
cone ind almost constant in most of the convergent part of the nozzle; at the 
e (Cc) throat it is still insignificant, and it could be compensated by a reduction 
b the of only 0-00416 in. in the core radius b, = 0-278 in. Thus the flow in the 


1pies main stream is almost unaffected by the boundary layer. 
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This work forms part of an investigation into swirling motion supported 
by the Mechanical Engineering Research Organization of the Department 
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NUMERICAL EVALUATION OF INTEGRALS OF THE 
FORM I= [ fixje# dx AND THE TABULATION OF 
x 
THE FUNCTION Gi(z) = (1/7) | sin(uz-+4u*) du 
ri) 
By R. 8. SCORER (Imperial College, London) 
_ Reec:ved 31 August 1949] 
SUMMARY 
s type of integral often occurs in the calculation of the wave form due to a 
in a dispersive medium. The principle of a stationary phase is used to evaluate 
rms of the Airy Integral, of which tables are published, and of 
I 
(al = sint(wz hus) du. 
\ table is given of this function fo1 0 and of a closely associated function for 
0. In a separate note by Dr. J. C. P. Miller it is shown that these, together 
th the published tables of the Airy Integral, provide the solution of the differential 
/ P(z) 
re P(z) isa] nomial or a power series. 


1. Introduction 
INTEGRALS of the form 
r= [ fiw)eis@ de 


cur in calculating the wave form in a dispersive medium due to a source 
vhich can be represented by a Fourier integral. They may often be 
evaluated by making use of the well-known principle of stationary phase. 
In this treatment such a method is assumed to be applicable. It depends 
ipon finding a, to satisfy 

dy = ¢'(ao) = 0, (1.1) 
vhere dashes denote differentiation with respect to x. If f(x) is slowly 
arying near a ry and ¢ is sufficiently oscillatory elsewhere, then 
nly in the neighbourhood of x, is significant contribution made to the 
ntegral. If (1.1) has more than one root in the range of integration the 
ontributions must be added, but we may suppose without loss of generality 
that there is now only one. 


The standard procedure is to W rite 


/ 1 (2) ( exp{i(dy ; (x X9)*h5)} dx, (1.2) 


(Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950) ] 
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and pr‘ 


which is expressed in terms of 
i near to 


| exp(ix") dx, 

whose value is known. Such an approximation is only valid when the where 
succeeding terms in the Taylor expansion of 4(a), viz. ' Here 
1 (x —aty)9h?” +... (2.1), ( 
are negligible over a sufficient range of x. The criterion as given by (3) anc 
Lamb (1) is that ” ” Ja 
ipo << po i. (1.3) 3, Ca 
The following method was devised for a special problem in which ie 

(1.3) did not hold. 


2. Formula for the integral 

When 45 = 0, a method using the Airy Integral given by Jeffreys (2) which 
may be used. That method is now generalized to apply when it is per- 
missible to write 


f(x) do T (x Xy)*ho T 4(x 5 X)*ho (2.1) and t 
- . now 
a+ba+ca?+dz 3, say, . - 
Sir 
a, b, c, d being thus known functions of x). Let purp 
| Pees ee a 
K(z) = Ai(z)+7 Gi(z) — | expi{a(uz+ 3u*)} du, (2.2) 
0 
then if two new numbers / and m are introduced, where 
since 
u l(a—m), ' 
inte: 
the exponent becomes a cubic in x which may be identified, apart from a 
constant, with i¢(7) as expressed in (2.1). For (2.2) becomes 
, lf {_ 1]3y3 3,2 ; m 3s ‘ To: 
K(z) = — exp{t| —2lm — 3 m3 + (2l+-Bm?)a—Bma?+4B25}} dx, (2.3) 
m 
and we are to assume that : 
O 


To 


I = f(z,) [ expfi(a +-ba+-cx*+-dx3)} da; 


1 


wherefore if also 


b 2l- jak i.e. 1 (3d)! | 
c Bm m —c/3d (2.4) 
d 4/8 | Zz (3d)-4(b- oe, na) 
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snd provided that either m is close to x, when 2p is near 2, or Xp is not 


near to vy Or Yo, WE have 


I 7 J (%o) iBK(z), (2.5) 

he ce 
‘ | ere 6 a+-zlm L/3978 a c— + - 2.6 
nth ; = 7 3d ° 9d? ( 
Here 7, 8, and z are functions of 2 that are readily calculated from 
9.1), (2.4), and (2.6) if d(x) can be differentiated. Ai(z) is already tabulated 

on by | (3) and so J is given by (2.5) when Gi(z) is tabulated. 


(13) 3. Calculation of Gi(z) 


The problem is thus reduced to that of evaluating 


whicl 
13/4 ( ; > 14,3 ¢ 
Gi(z) — | sin(uz+3u ) du, (3.1) 
vs (2 hich is a particular integral of 
; Der ai ] 
7 ay ee (3.2) 
(2.1 nd this equation may be integrated numerically if Gi(0) and Gi‘(0) are 
Known. 
Since Gi(z) oscillates for z 0 it is more convenient for interpolation 
purposes to tabulate, for z 0, the function 
(2.9 Hi(z) lf er ‘ 
2.2 1i(z) — | exp(uz— ju ) du (3.3) 
Bi(z) —Gi(z), (3.4) 
nce Bi(z) is already tabulated in (3). It is readily shown by contour 
tegration using the method indicated by Miller (3) that 
om a : ; 
Gi(0 + Hi(0) LBi(0) 3-*Ai(0) 3-5/(—4)! (3.5) 
Gi’ (0 1Hi'(0) = 4Bi'(0) 3-1Ai’(0) = 3-#/(—2)! (3.6) 
— To twelve decimals the values are 


Gil0) 0-20497 55424 78 (3.7) 


Gi'(0) = 0-14942 94524 49. (3.8) 
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was used to determine the values of Gi(z) for the first two steps. Noy 
(3.9) can be written 


, : ij . & . 3.6¢ 
Gi(z) = « Ai(z)+ 8 Bi(z) —— |— + — +—-+ ---]» 
where « and f are numerical constants. 
Hence 


Gi(0) x Ai(0)+-B Bi(0) 
Gi'(0) = a Ai’(0)+-f Bi'‘(0), 


and comparing these with (3.5) and (3.6) we find « = 0, 8 = 4, Le. 


a2 925 ‘ 28 
ie 32 pe (3.10) 
2 ae 


The integration was performed by using the formula 


: mn | 
2 1, 31 ps... FY (3.11 


é°y (2) 14 o- — o* + — c =e 

12 240 60480 dz* 

The values of the differences on the right-hand side were obtained at each 

stage where necessary by performing the integration a few steps ahead as 

necessary with progressively fewer significant figures and thus also fewer 
terms in (3.11). 

The functions were tabulated in the range —10 < z < 10, outside which 
the asymptotic formulae give at least 9 correct decimals. Ten decimals 
were used in the integration and the asymptotic formula was used at 
Z 8, 9, 10 together with the series (3.9) at z 5 to eliminate rounding 
off errors which grow like Bi(z) in forward integration for z > 2. Forz < 0 
the unwanted solutions of the differential equation (4.2), below, satisfied 
by Hi(z), are oscillatory and do not grow unduly. A programme had been 
made for tabulating Ai(—z), by integrating the differential equation which 
it satisfies, on the EDSAC at Cambridge; and no major modification was 
required to tabulate Hi(—z). The table of Hi(—z) was therefore made on 
the EDSAC using an integration interval of 0-02, every fifth value being 
printed. The method has been described by Wilkes in (4). 

To ensure that the oscillatory errors were small the integration was 
carried to z 12 and compared with the results obtained by hand for 

10 >2 12 starting from the asymptotic formula. Eight decimals 
correct were obtained throughout, and a larger table is contemplated to 
contain these; but for present purposes, where only modified second 


differences are given for interpolation, seven decimals are set down. 
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No 4, Asymptotic formulae 
The asymptotic form of Gi(z) can contain no multiple of Bi(z) when 
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- © since Bi(z) increases indefinitely. Now a particular integral of (3.2) is 


Lyi .-2!.. 8} 8! 
a 2A" 327°" 3.6210 TOP 
hich is of a higher order of magnitude than Ai(z) whose asymptotic 
rm is ‘ | 
=—— 2-* exp(— §2°)(1—g2-*+-...). 
a NTT 
Thus, in the Poincaré sense, for z 0 
” l/l 2! 5! 8! 
Gi(z) ~ cs . oe * (4.1) 
7 a" 327 © 3.6210 | 
(3.10 ; : er . 
and the numerical work has confirmed that no multiple of Ai(z) enters into 
the formula. 
When z 0, writing z we have (see (2) or (3)) 
9 ] ° . 4 ° « 
(3.11 Hi(¢) = Gi(¢)— Bi(¢) exp(—ul—4u?) du. 
0 
hes Integrating by parts, the first two or three terms of the asymptotic 
“ad as expansion of the right-hand side can be obtained, the subsequent terms 
lew being found by substitution in the differential equation satisfied by Hi(z), 
umely > 
whic : zy (4.2) 
ima ; ° F = 
e Thus, for 
. — - ijl 62! 6S 8! 
din Hi Bi(f)—Gi(Z) ~ : — = a re P (4.3) 
( 7TA\C ( 3G! 3.6220 
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NOTES ON THE SOLUTION OF THE EQUATION 
y'—xy = f(x) 


By J. C. P. MILLER and ZAKI MURSI 
(Scientific Computing Service, 23 Bedford Square , London, 
and Farouk I University, Alexandria, Egypt) 


Received 31 August 1949] 


SUMMARY 


By suitable approximations and changes of variable, many second-order linear 
ferential equations may be reduced to the Airy Integral form 
| a 7 | 
ry J (2). 
l 
hese notes it is shown how this equation may be solved when 
j du 
1 (x) Oxr).u a 
u a Ai(x) +b Bi(x), 
1 where g(a) and h(a) are expressed as power series. The solution may be expressed 
same form, or as a series of derivatives of u. The solution is also given in the 
e where f(x) is itself expressed as a power series; in this ease it is of the form 
y h(a U(x).v, 
u mor) 
which | and 1 ire expressed as power series. The function Gi(a) is tabulated 
previous paper (Scorer, 1). The series in the solution terminates if g(a) and h(a) 
polynomials, or if f(a) is a polynomial. 


l. Ir is well known that many problems in mathematics and physics 
depend on the solution of linear differential equations of the second order. 
112  Byappropriate choice of dependent or independent variable, or both, such 
equations may be reduced to the normal form, that is, to a form not 
containing a term involving the first derivative of the dependent variable. 
We may therefore, without loss of generality, write the equation to be 


solved in the form dy 


+. I(x).y F(x). (1.1) 
3 dx” 
Over an appropriate range of the independent variable we may in 
certain cases expand the coefficient /(x) in a power series 


I (2) I f, (a Xo) +-I, . (w—2Xp)?/2!-4 eeey (1.2) 


and an approximation to the solution, sufficient for many purposes, may 


yt 


ve found by neglecting terms other than the first one, two, or three. It 
may be possible to divide the full range of x for which a solution is desired 


[Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 


9092.9 


I 














































114 





J. C. P. MILLER AND ZAKT MURSI 


into a number of intervals, each such that a one-, two-, or three-term B 
expansion (1.2) is an adequate approximation. The solutions valid within | 


the various intervals may then be fitted together at the ends of the | whe 
intervals, and an approximation obtained covering the whole range of x ° 
under consideration. ‘ 

If we take only the first term of (1.2), the function J(2) is replaced by whe 
a constant in each interval of 2, that is, by a step-function. In each oe 
interval, (1.1) becomes the familiar equation with sine and cosine, or with 
positive and negative exponentials, as complementary solutions. | 

If we take intervals in x such that J(x) may be interpolated with second 
differences only, third differences being negligible, then a quadratic 
approximation to /(x) is taken, and the equation (1.1) has a comple- 3. 
mentary solution in terms of Weber functions. 

This paper is concerned with cases where intervals of x are chosen so its 
that linear interpolation of J(x) is sufficient, i.e. so that only the terms ” y 
I, + [,(2—ay) are retained in (1.2). In this case the complementary function |“ 
in (1.1) is expressed in terms of the Airy Integral Ai(é) and the associated 
function Bi(é), by writing then 

L,+(fh—1,%) = —Hfé. as) 7 
The equation (1.1) is theh converted to the form ” 
: * whe 

wat f(®). (1.4) 

This is the equation dealt with below, though x will be used in place of . 
then 

2. Using suffixes to denote differentiation with respect to x, the equation ite 
to be solved, for certain forms of f(x) is ; 

u 
Yo—xy = f(x). (2.1) oe 


This suffix notation for differentiation will be applied only to the variables —x*u 
y and u below. te 


Before dealing with the equation (2.1) it is convenient to develop certain =u : 
properties of the complementary function. We write , Pu 
d? TI 
od —s = eae Ay (2.2) 
dx 
for the operator on the left of (2.1), which becomes dy = 0, with the 
general solution . 9 § 
F u = a Ai(x)+5b Bi(z), (2.3) 
in which a and 6 are arbitrary constants, while Ai(x) and Bi(x) are the vv 
; . P pe : as ae givin 
Airy Integral functions, as defined by Jeffreys (2, pp. 476 and 477) and like | 
: “ 


tabulated in (3). 





cond 
lratic 


nple 


en. so 


> 
erms 


ctiol 


iated 


(1.4 


» ol é 


1atiol 


(2.1 


‘jables 


ertall 





By differentiation it is readily verified that 


u > TU 


+nU,»-» (nm > 0), 
whence w, can be expressed in the form 
u, DP, (@)-U+G, (2X). U4, 


where Py (%) and q,\%) 


Thus, 
Us xU Us (x? +-4)u-+-62ru, 
! Oy2 “3 _| 
Us u a Uy Un Qxr-u + (a = 10)u, 
Us x*ut2u, Ue (x4+-28x)u+- 12274, 
U- fru LU, Ug (1623-4 28)u- (a4 52x)u, 


3. The relation (2.4) may be written as 


vu nu, 1 (n - 0), 
which gives an immediate solution to (2.1) when f(x) 
U,,/N. Thus if 
(x) Ay UA, Uy + Ag Ug Az Ug 
then y Ay Uy 5a, Us- 4A Us ' las i eee 
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are polynomials of degree not exceeding }n+1. 


- 


(3.1) 


U,-1, namely 


(3.2) 
3) 


Co 


The form (3.2) for f(x) is less unlikely than at first it appears to be. 


For if f(a) 


may be expressed as power series, 


g(x).u+h(a).uy,, 
vhere q(x) and A(a) 
JX) = Jo 
h(a) hy +h,x- 
then, by application of (2.4), 


f 


of x, we find 


(3.4) 


(3.5) 


repeated if necessary, to eliminate powers 





u Uu Uy Uy 
i] Us LU, Us Uu 
l u 2u r-u U-—4U. aise 
4 l 1 5 2 (3.6) 
1 U,— 6U,-+ 2u vu; Uy 9U, + 8u, 
l= Ug—12u,+ 20u vu, U,— l6u,+ 44u,—8u 
‘ o> a 
l U;9— 20u,+ 80u,—40u, CPU, U4, — 29Ug 4 140u,; 140u, 
Thus (3.4) and (3.5) may be expressed in the form (3.2), where 
Ay Jot 293+ 409,+...—h,—8h,—280h,—... | 
») » 
a, 29.— 409; hy +8h,+ 280h,+-... (3.7) 
Ao q; 209, th, 140h, tue. 
Ay 6g,—...+4,+-44h,+-... 
giving the solution (3.3), which may, if desired, be converted to a form 


Ke (3.4) by use ol (2.6). 
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This method of solution is most convenient when g(x) and h(x) are 
polynomials, or may be treated as such. 


4. Next let us suppose that 
f(x) b,-4 b,x | by a* } bg a+... . (4.1) 


To solve (2.1) in this case we note that 


tan n(n—1)ar—-2#—gnrt 
and, in particular, that 
‘a x 
ax a (4.2) 
ox xt2 
ta -x44+-6x 


while it is readily verified that 


v a Gi(x)+u, (4.3) 
- Rr .. — 
where Gi(x) = sin(wa+4u) du, (4.4) 
“3 
satisfies the relation ov ib. (4.5) 


The equations (4.2) may be rewritten to give single powers on the right, 


thus 
H(—v) e 3 
#(—1) + | 
o(—x) = 2? (4.6) 
o(—x?—2v) = x3 | 
o(—a?—6) x4 


and so on, the general expression being 


H{ —x"—n(n—1)a”-3§ —n(n—1)(n—3)(n—4)a”-6—.,,,} gett, (4.7) 


The series in the brackets terminates if n+1 is of the form 3k-+1; if 


n+1 = 3k, the series continues to the term in x, and is followed by a 
final term 
—~n(n—1)(n—3)(n—4)(n—6)(n—7)...5.4.2. 1.0. (4.8) 


The solution to (2.1) with f(a) given by (4.1) is thus 


u—v(bo+1.2b,+1.2.4.5b,4+...) 
—(b,+2.3b,+2.3.5.6b,+...) 
-a(bo+3.46,+3.4.6.7b,+...)— 
—2*(b,+4.5b,+-4.5.7.8b,+...)—... 


y u—b,v—b,—b,x—b,(u?+- 2v) 4 


(4.9) 
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are 5. The solution (4.9) may be written as 


Y u Cv dy, d, 


x—d,x*—d,x°—.... (5.1) 
Then, substituting in (2.1) and using (4.1) and (4.2), and equating coeffi- 
1] cients of powers of 2, 

2d,—by—c = 0) (5.2) 
(n+1)(n+2)d,,,.—b,—d,_, = 0)" zi 
These equations do not suffice for the determination of all the coefficients 
but if c, dy, and d, have been found, e.g. from their explicit forms in (4.9), 

1.2) " then (5.2) may be used to obtain further coefficients in succession. 
Again if f(a) is a polynomial of degree k we may assume d,, 0 for 

k+-2, and determine coefficients d;.,5, d;.,, d;,, etc., in succession. 


6. More generally, by writing (2.1) in the form 


a (x Yo : 

13) » y — 1") _¥ (6.1) 
: x x 

‘44 nd using successive approximations where appropriate, we obtain the 


formal solution 


én . Ld fxe) laf & =) 1 d* {1 ral; d* 4} 


x" vr dx x a da*\x da* x x dx? la da*\xa du x J 
iht, (6.2) 
Similarly we may write Jo = f(x)+2y, (6.3) 


vhence the formal particular integral (in a condensed notation) 


(4.6 y= || f(x) (dx)? || x || f(x) (de)t+ |] aw |] x || f(x) (dv)o+.... 
(6.4) 
It is interesting to note that, if f(x) 1/7, equations (6.2) and (6.4) give 
1/1 1.2 1.2.4.5 (6.5) 
~ Y~ t = " re do 
(4.7) : = yi a 
i au 1 / a? Yr a bel 
and Yy | = ee (6.6) 
by a 7w\l1.2 1.2.4.5 L248. 4.8 
' hich are the expansions fol Hi(a) Gi(x) 31 (2) and Gi(x) —1 Bi(x) 
.o) . 
\ respectively (see (1) (4.3) and (3.8)). 
7. These ideas have been much more fully worked out in an unpublished 
Ph.D. thesis (Mursi (4)), where extensive tables of coefficients and further 
related developments are given. This thesis is inaccessible, having been 
Q) } - 7 ° - > ° . 
(4.: submitted to Farouk I University, Alexandria, Egypt. Various properties 


f the function Ai(2) {4 Bi(a)—Gi(a)! are also discussed; several results 


en in (1) have been obtained independently. 
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INDUCTION OF ELECTRIC CURRENTS IN 
UNIFORM CIRCULAR DISK 


pp. | 
By A. A. ASHOUR 
tae (Math. Dept., Faculty of Science, Cairo University) 
Received 13 July 1949] 
; SUMMARY 
1 by By regarding a uniform disk, or any symmetrically conducting surface of revolu- 
465 * tion, as composed of an infinite number of co-axial annular circuits, the determination 
f the electric currents induced by an external field is reduced to the solution of a 
Fredholm integral equation whose kernel K(a, x’) has infinities of order log(a—.’) 
en a v’. Two methods of solving this equation are described and illustrated. 


time constant for a uniform disk is calculated and agrees approximately with 


value estimated by Lamb. 


» 1, Introduction 
THE induction of electric currents in a conducting circular disk by a 
ying magnetic field was discussed by Lamb (1,2) by treating the disk 
is the limit of a thin spheroidal shell when its axis becomes vanishingly 
small. In this way he obtained solutions for those cases in which the 
resistance at distance r from the centre of a disk of radius a is a special 
function of (a2—r?)!, but he was unable to obtain the solution for the case 
when the resistance is uniform. He was able, however, to estimate a lower 
limit for the principal time constant of a uniform disk by regarding it as 
in intermediate case of two known solutions. As far as the writer knows 
there has been no further mathematical treatment of the uniform disk 
problem, though Bruckshaw (3) has made an experimental investigation 
of the resulting induced field, with the object of obtaining results useful 
1! geophy sical prospecting. 

In the present paper a method suggested to the writer by Professor 
\. T. Price is developed for treating the problem when the conductivity 
of the disk and the inducing field have axial symmetry. The disk is 
regarded as composed of a large number of concentric annular circuits, 
and the problem thereby reduced to the solution of a Fredholm integral 
equation. Two methods of solving this equation are described and illus- 
trated with a numerical example. The time constant for a uniform disk 
s calculated and agrees approximately with the value estimated by Lamb. 


2. Derivation of the general equations 

A thin circular disk of radius a and uniform surface resistance p per 
init length is situated in a varying uniform magnetic field H. If the field 
8 oblique to the disk, only the normal component will be effective in 


(Quart. Journ. Mech. and Applied Math., Vol. III, Pt. 1 (1950)] 
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inducing currents in it, and these will be symmetrical about the axis. It | Hen 
is therefore sufficient to treat H as normal to the disk. 

The disk will be regarded as composed of a large number of concentric 
annular circuits, each of inner radius 7 and outer radius r+-dr. The total 


, . 
' ; . “i 1i¢ 
resistance of the ring (r,r+é6r) will be 27rp/d5r and the total current * 
circulating in it is /(r) dr, where J(r) is the integrated current density, 
If N(r,t) is the total induction through the circuit at time ¢, we have ; 
Now 
C y 
2arpI(r) —— IV (P, $). (1) 
ot whet 
If M(r,r’) is the coefficient of mutual induction between the two rings r his 
and r’, N(r,t) is given by are ¢ 
Nir, t) nr?H 4 > I(r’) dr’ M(r,r’). (2) The 
‘i ; 
The summation in (2) covers all the circuits including r’ =r. If dr is —— Jog( 
infinitesimal, (1) takes the form: > then 
a Vv ; anni 
r ) , 217, Tf , ‘ 
Ir) = —2 pH [or ar’, (3 
2p 2rp . 7 
: Thu 
where p denotes the operator (é/cet). If p is treated as a parameter, equa- whe 
tion (3) is an integral equation of Fredholm type to determine /(r, p). 
4, 7 
3. Extension to surfaces of revolution and to axially symmetric T 
inducing fields cent 
This method can be applied to the more general case of a thin non- so t] 
uniformly conducting surface of revolution, whose conductivity has axial the 
symmetry, situated in a varying magnetic field which is also symmetrical 
about the axis. Let the equation of the surface be 
z= f(r), (4) whe 
and let the z-component of the inducing field be F(z,7r). Consider the Ea 


annular circuit with radius 7 and breadth ds. The resistance of this circuit 
is 27rp(r)/ds, and the current flowing in it is J(r) ds. The part of NV due nan 


to the induced currents will now be given by 


si r whe 

( Mi(r,r’)L(r’) ds’ | M(r,r')1(r' A+f 2(r)}t dr’, (5) 

1 r ry 
where 7,, 7, are the radii of the boundaries of the surface. The part due Che 
to the inducing field will be “4 

: 6 
[ 2ar’ F(z,r’) dr’ = mr*H(r,t), say. (6) Ry 


0 
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s. It Hence i 
pH | *, t) , | , ’ 1 "79 , 7 — 
I(r) dtd E I(r’) M(r, rf (rH ar’, (7) 
ntric 2p(r) 2rrp(r) P : ‘ 
r1 
total 


which is a Fredholm integral equation with kernel 


rrent 
sit 1+ f(r’) PM (r, r’)/2arp(r). (8) 
| , |: oo eo 2 yt 
Now M(r.r’) t7r,/ (77 ) ( e) K(e)— E(c)}, (9) 
| \e ony 
| vhere c drr’S(r+-r’)?+ h?\-1, 
ngs h is the normal distance between the planes of the circuits, and K(c), E(c) 
wre complete elliptic integrals of the first and second kinds respectively. 
2 The coefficient of mutual induction will be finite everywhere except when 
1. This is the case when 7 vr’ h 0. In this case H(c) behaves like 
or Is yo(r—r’) and M(r,7r’) reduces to the coefficient of self-induction ; we can 
* then use the well-known approximate value for the self-induction of an 
nnular ring of inner and outer radii 7, r’ respectively, namely 
5) M(r,r’) t7r7| log{8r/(r yr’) il. (10) 
Thus the kernel of the integral equation (7) is finite except atr = r’,h = 0, 
equa here it becomes infinite like log(r—r’). 
4. The uniform disk 
etric The solution of the integral equation (3) will now be considered; the 
general equation (7) can be treated on parallel lines. If we write x = r/a, 
non so that M(r,r’) = aM(x,2’), and also write J(x) for I(r), equation (3) takes 
axia the form : 
— I(x) = wH +“. | 1(x’)M(w,2') de’, (11) 
ae 
t vhere A ap/2p. (12) 
r the Equation (11) can be transformed to another with symmetrical kernel, 
‘ircult 1 
V du namely us(a) AHa tA | ob(x’)K (a, x’) da’, (13) 
0 
ee us(a) t I(x), 
5 (14) 
K (x,a M (x, x") /{47r4/(xx’)}. 
rt due The new kernel K(x, 2’) will still have an infinity of order log(a—a’) when 
Uw « 
One method of solving (13) approximately is to regard the disk as 
(6 composed of, say, ten circuits each of breadth (a/10) and with mean radii 


V05a, O-l5a, 0-25a,..., and 0-95a respectively. When calculating the 
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coefficient of mutual induction between any two of these circuits they are 
treated as line circuits. The expression (10) was not, however, used to 
calculate the coefficient of self-induction for the above circuits, because 
the breadth of each circuit (especially those nearest the centre of the disk) 
although small, cannot be neglected in comparison with the mean radius 
of the circuit. Instead each ring of mean radius 7 was further subdivided 
into ten annular rings, each of breadth 0-0la, and of mean radii r,, r, 


and 719, where r, = r+(0-01s—0-055)a. (15) 


If M(r,r,), M(r,79),..., M(r,r,) are the coefficients of mutual induction 
between the ten new rings and the original ring, its coefficient of self- 
induction is approximately 


6 2, M(r, ¥.}. (16) 


n 


Using the values for log{M(a,x’)/47,/(vx’)! given by Maxwell (4), the 
values of K(x,x’) are readily obtained from (14) and (16). 


5. The solution of the integral equation by the method of ‘con- 
tinued substitution’ 

In this method, which is a special case of a general method due to 
Price (5) for any non-uniformly conducting sheet, a first approximation 
ws, (x) is obtained for ys(a) by neglecting the second term on the right-hand 
side of (13). Thus we have 

w(x) = AH! = AH¢,(z). (17) 
This approximation corresponds to neglecting the self-induction of the 
disk. The second approximation is obtained by substituting y%,(2x’) for p(x’) 


in (13). We get yb5(2) AH¢,(x) | 2Hd,(2), (18) 
1 

where d(x) = 4 | by(x')K (x, x’) dx’. (19) 
0 


By continuing this process, the nth approximation will be given by 


wp, (ac) = AH > _y(x)a™-2, (20) 
m=1 

; 

where bie) = 4 | bu —1(x' K(x, x’) dx’. (21) 
0 

Hille and Tamarkin (6) have shown that the equation 
b 
p(x) = f(x)+A | K(x, x’)b(2x’) da’, (22) 


a 


where K(a,2x’) is symmetric and has discontinuities when x = 2’, admits 








1 con 
that 


Hene 
to inf 


It fol 


Th 
0-95 1 


above 


This 
repre 
used 
operé 
state 
can t 

As 
diseu 
of ra 
em. 


moni 


A 


y are 
» to 
ause 
disk 
dius 


vided 


atio1 


-hand 


(17 


yf the 


/ 
aie a 


(18 


(19 


(20 


(21) 


» 


dmits 














ELECTRIC CURRENTS IN A UNIFORM CIRCULAR DISK 123 





convergent solution by the method of continued substitution provided 


nat 


(i) f(x) is integrable (in the sense of Lebesgue) ; 


} 


l nes ; 
(ii) K (a, x’) dx : . 
b a A 


jence, since x! is integrable, the infinite series obtained by making n tend 
infinity in (20) will represent the solution of (13) provided that 


1 


1 

A 4 K(x, 2") de| (23) 

t follows that I(x) = d(a)a AHx-? ¥ ¢,,(x)A™. (24) 
m=0 

The functions ¢,,(x) are evaluated at the points x = 0-05, 0-15...., and 


95 using their recurrence formula (21), and the values 2~!¢,,(”) for the 


bove values of 2 and n 1, 2...., and 8, are tabulated in Table I. 


TaBLE [. Values of x 








15 55 0°65 O75 0°55 o°95 
5 4 0°65000 | 0*75000} 0*55000)| 0*95000 
I 2 5641 2°2199 375078 3°4889 3°1616 
I 767 15°517 15°679 14°77! 12°647 
7 70°075 | 69°955 | 94°492 | 54°340 
5°405 I 29 75°55 305°06 318-33 312°53 285°83 | 239°67 
7 1385'8 1435°5 1402°5 1288°0 1069°5 
737°3! 3 i! } 37°5. | 6261-4 | 5465°9 | 6306°5 §737°3 | 4797°7 
25235 29143 25304 25503 21570 
17 127260 | 131240 | 127800 | 116150 | 97077 





This table actually contains the coefficients of the powers of A in the series 
representing the solution of the integral equation, and hence it can be 
used in solving any particular problem with any special value of the 
operator A. In case of periodic external fields, if we only require the steady 
state solution, we may write H as the real part of He’™; the operator 0/é¢t 
can then be replaced by 10. 

\s an illustration, a numerical case is considered which is of interest in 
discussing the effects of a large area of sea water in geomagnetism. A disk 
i radius 1,000 km., thickness 5 km., and specific conductivity 4 10-" 


ay 


n.u. is situated in 


a uniform field normal to it and having simple har- 
monic variations of period one day, say H 


0-07268:2 


H cosat. These values make 


The amplitude and phase of the current density, which 
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is represented by the real part of the right-hand side of (24) are cale 


A. A. ASHOUR 


ulated This st 













and tabulated in Table II. K(x, x 
TABLE II functic 
= rs an VIA). \ 
x 0°05 O15 0°25 0°35 O45 0°55 O05 | O75 0°55 : 
R 4 | oror102 | o-or7¢ 0238 | 0:0308 | 0:0377 | 0-447 | o°0519 | o-o5g1 agon 
x(x) 65°8 Sy 65°5 66°4 67 °c 67°9 69°4 | 711 73° D anc 
The amplitude and phase of the induced current density tively. 
I(x) HR(x)cos{ at x(a)} e€.m.u. discon 
The above table shows that the amplitude increases steadily as differ : 
the distance from the centre increases, till it attains its maximum at the’ “ Le. 
boundary. On the other hand, there is little change in the phase of the 
induced current density. If the self-inductance of the disk be neglected, 
the induced current density will be given simply by 
I(x) vA, (25) | Henee 
functi 
so that R(aw) = x and a(x) = 7/2. hee 
teristi 
6. Fredholm solution of the integral equation Hill 
Fredholm (vide (7)) has shown that the unique solution of the equation _ ties o 
b that 
(a) = f(x) +A | ba’) K(x, x’) dx’, (26 
a 
when the kernel is continuous and finite for all x, 2’ in the interval (4,5), — In thi 
is given by b provic 
S Pee ; , a 
(a) = f(x) + | f(a’)V (a, 2’, A) dx’, (27) used. 
V(A) . The ¢ 
a 
where numb 
4 ( 1)” , * ( l n ‘ 
V(A) > —— a A V(x v’. A) S yb, (2 a’ )Ant 
i —_ nt 
n=0 n=0 ' 
b b As far 
a | : | D, dx, Ax... AL», SO. si] 
a a are ol 
b b F 
: ; " 
b, [2 3 | 2 D(x, x ) dx, Efe.cAEa, 
a a 
D,, denotes the determinant A{K (x,, 2,)}, (25) 
| 
r Re desea n; 8 Se Mee nN, 
and Wher 


K(x, x’) 
K (2,, x’) 





lated 


ut the 
Tt the 


Pcter 
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this solution is not immediately applicable to (13) because the kernel 
r,2’) has an infinity at 2 = 2x’. Bdocher (8) introduced the modified 


netions V(A), V(a,a’,A). These are obtained from the solving functions 


VA), V(w,a’,A) by putting K(x,,2,) = 0 (r = 1,2,...,m) in the principal 


wonal of D, and D(x,x’). The corresponding values for a,, 5,(x, 2’), 
) and D(a, 2’) will be denoted by @,, ,(x,x’), D,, and D,,(x, 2’) respec- 
ely. Bécher then proved that for any symmetrical kernel whose 
scontinuities are regularly distributed, the modified functions do not 
ffer from the solving functions except by a constant exponential factor 
1.e. 


V A e c)V() 9 
(A) xp(c)V(A) (29) 


V(x, x’, A) exp(c)V(a, x’, A). | 

Hence for such kernels the solution of (25) is not affected if the modified 

nections are used instead of the solving functions. Also, the charac- 
ristic numbers of the kernel are the zeroes of either V(A) or V(A). 

Hilbert (9) has extended this theory to symmetric kernels having infini- 

es on the line a x when a positive number B < 4 can be found such 


+ 
it 


lim (a—a’)PK (a, 2’) = 0. (30) 


In this case he proved that Fredholm’s solution will be valid and unique 
rovided that instead of the solving functions the modified functions are 
ised. With the help of the latter functions all infinities are removed. 
lhe condition (30) is satisfied by equation (13), since for any positive 


) 
imbet Db 
, 


lim (a x log (x x) Q. 


, 


\s far as the writer knows no simpler forms for @,,, 6, (x, x’) have been found. 
So, since this is necessary for any practical application, the following forms 
ire obtained for any symmetric kernel. 


From the definition we have by expanding D,, 


‘ . 1)r P _ . 
a (n mS! LG... #> 1, @ l, a, 0, (31) 


i. ( { K(x, to) dx, dx 


} 


» | A(x, Xy)K (aq, %5)...K(X,_1,%,) K(x,, 2) dx, dx,...dx,. 
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Similarly b, (a, 2’) = n! 3 — soe 4 oo (33) 
* n (n r)! r n=? 
where Iy(x, 2’) K (x, x’), 
b 
I (2, 2’) | L_s(z, €)A(é,2’)d& (r > 1). (34) 


The following two useful relations follow immediately : 
bb 
Lig = || K(x, 2%')E(x, 2x’) dada’, | 


aa 


bb (35) 

- -(r+1) | | K (x, 2’ )b,, (x, 2’) dada’ | 
aa 

All functions of two variables needed for our special case have been 

evaluated for the hundred points (x, 2’), where 2 = 0-05, 0°15.,..., 0-95 and 

x’ 0:05, 0°15,..., 0-95. The series V(A) and V(x,2’,A) are evaluated as 








far as the term A‘. Starting with J,(v,2’) = K(x,x’), the functions J, (z,2’), | 
n 1, 2, and 3 are evaluated by carrying the integration in (34). The 


integrals J,, J;, and J, are then evaluated using (35). Their values are 
found to be 
i, = 1-8161, I, = 1-6404, I, = 1-7018. (36 
The coefficients d,, were then evaluated, giving 
V(A) 1 —0-9080A2— 0-5468A3— 0-0132A4... . (37 
It will also be seen that 


b,(z, 2’) = I,(2, x’), b, - I, Ig(2, x’) + 21,(a, x’), 


ae,2°) = Lfs,2"), b, — 21, L(x, x’) +31, L(x, x’)—61,(x,2x’), (38 


and thus the first four terms in V(a,2’,A) can be obtained readily. Thi 


solution of (11) can now be obtained in the form 


1 
4\NH ff ooo 
I(x AHa +. —. v'iV(a.a’.d) da’. (39 
(x) x!V(A) V(a r ) da 


0 

The special numerical case already discussed in section 5 was als 
calculated by the present method. The amplitude and phase of th 
current density are obtained at the same points as before. The results 


given in Table III, agree well with those in Table II. 
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93) ras.e IIT 
"45 0°55 0°65 o°75 0°85 0°95 
34 ‘0170 "0308 | 0°0377 | 0°0447 | O°051g | O'0591 | 0°0664 
6a°5 Of 66°9 67°9 69°4 7I°! 735 76°6 





e amplitude and phase of the induced current density calculated by Fredholm’s method. 


7. The free currents in the disk 

If a symmetrical set of currents is excited in the disk and left to decay, 
there being no inducing field, equation (13) will be homogeneous. Assum- 
ng the current density everywhere to be decaying exponentially like e—”, 
the operator @/ét may be replaced by —7~! and A is then equal to 2a/rp. 
[he characteristic numbers of the kernel thus give the values of the time 
mstant rT. 

Lamb’s discussion of this problem led to the fixing of a lower limit to 


e been, the principal time constant, namely + = 2-26a/p. This means an upper 


95 ar mit to A = 0-885. This value should correspond to the smallest zero 
ted as} of V(A) which is given approximately by (37). Using the five terms of 
t PI y 2) g 

(x.a’),| V(A) given in (37), we obtain A = 0-850, which gives +t = 2-34a/p. This 


). The value, which is probably in error in the last figure, is slightly larger than 

ies are Lamb’s estimate which was a lower limit; the agreement is satisfactory 
und affords some check on the calculations.+ 

(34 It is of interest in geomagnetism to estimate the rate at which free 
currents decay in the oceans. If an ocean is regarded as a large uniform 


disk with radius a and depth d, the time constant, i.e. the time required 


(97). for the current density to reduce to a fraction 1/e of its initial value, can 
be calculated by using the above value of 7. This is calculated for the 
values of a and d shown in Table IV, taking the conductivity of sea water 
to be K $< 10-1! e.m.u. 

TABLE IV 

(38 

3 4 5 
T] ; 
2°34 
I 3°12 4°16 
(39 I )0 5°20 6°50 
Che time constant (in hours) for a circular ocean 
of radius a 1000c km. and depth d km. 


In conclusion | wish to thank Professor A. ; 2 Price of the Imperial 
alk College, London, for suggesting this investigation and for his interest and 

esults tee 

help in its development. 


One of Lamb’s actual problems will be discussed by this method in a following paper. 
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